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Abstract. We review the application of field-theoretic renormalization group (RG) 
methods to the study of fluctuations in reaction-diffusion problems. We first investigate 
the physical origin of universality in these systems, before comparing RG methods to 
other available analytic techniques, including exact solutions and Smoluchowski-type 
approximations. Starting from the microscopic reaction-diffusion master equation, 
we then pedagogically detail the mapping to a field theory for the single-species 
reaction kA — > IA (£ < k) . We employ this particularly simple but non-trivial system 
to introduce the field-theoretic RG tools, including the diagrammatic perturbation 
expansion, renormalization, and Callan-Symanzik RG flow equation. We demonstrate 
how these techniques permit the calculation of universal quantities such as density 
decay exponents and amplitudes via perturbative e = d c — d expansions with respect 
to the upper critical dimension d c . With these basics established, we then provide 
an overview of more sophisticated applications to multiple species reactions, disorder 
effects, Levy flights, persistence problems, and the influence of spatial boundaries. We 
also analyze field-theoretic approaches to nonequilibrium phase transitions separating 
active from absorbing states. We focus particularly on the generic directed percolation 
universality class, as well as on the most prominent exception to this class: even- 
offspring branching and annihilating random walks. Finally, we summarize the state 
of the field and present our perspective on outstanding problems for the future. 
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1. Introduction 

Fluctuations and correlations in statistical systems are well known to become large in 
the vicinity of a critical point. In this region, fluctuations have a profound influence on 
the macroscopic properties of the system, leading to singular thermodynamic behavior 
characterized by universal critical exponents and scaling functions. These power- 
law singularities can be traced to an underlying emerging symmetry, namely scale 
invariance: at the critical point, the system possesses a diverging correlation length. 
Therefore renormalization group (RG) methods, which explicitly address the behavior of 
physical observables under scale transformations, have been employed with considerable 
success in describing critical fluctuations. The renormalization group provides a natural 
conceptual framework for explaining the occurrence of critical behavior, the emergence 
of universality, and the classification of different systems in terms of universality 
classes. Moreover, RG tools (especially in conjunction with series resummations or 
numerical implementations) also enable quantitative, controlled calculations of universal 
properties. 

Most successful applications of the renormalization group address systems in 
thermal equilibrium, where the Boltzmann-Gibbs probability distribution provides a 
solid foundation for explicit calculations. However, many systems in nature cannot 
be cast into even an approximative equilibrium description, and a large variety of 
nonequilibrium systems, both relaxational and driven, also exhibit critical behavior, 
which should presumably also be able to be analyzed by RG techniques. Unfortunately, 
even for nonequilibrium steady states, we presently lack a general statistical framework 
to construct the corresponding probability distributions and hence obtain the relevant 
macroscopic quantities. Consequently, there are relatively few cases where such explicit 
calculations can be developed, at least to date. For systems that can be represented 
in terms of stochastic partial differential equations of the Langevin type, there exists a 
well-established mapping to a field-theoretic representation pQ. However, this inherently 
coarse-grained, mesoscopic approach relies on an a priori identification of the relevant 
slow degrees of freedom. Moreover, far from equilibrium there are no Einstein relations 
that constrain the generalized stochastic forces or noise terms in these Langevin 
equations. Thus, although the functional form of the noise correlations may crucially 
affect the long-time, long-wavelength properties, it often needs to be inferred from 
phenomenological considerations, or simply guessed. 

Fortunately, in certain cases an alternative approach exists which allows these 
fundamental difficulties to be at least partially overcome. This method relies on the 
introduction of a 'second quantized' ladder operator formulation IS] for certain 
classical master equations, and on the coherent state representation to construct the 
statistical path integral j3j. This in principle permits a mapping to a field theory 
starting directly from a microscopically defined stochastic process without invoking any 
further assumptions or approximations beyond taking the appropriate continuum limit. 
In this way, the difficulties of identifying the slow variables, and of guessing the noise 
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correlations, are circumvented, and a field theory can be straightforwardly derived. 
Subsequently, the entire standard field-theoretic machinery can then be brought to bear, 
and progress made on understanding the role of fluctuations, and on identifiying and 
computing universal quantities. (As a cautionary note we add, though, that the above- 
mentioned continuum limit may not always be trivial and benign, and occasionally 
additional physical insight needs to be invoked to obtain an appropriate effective 
field theory description.) In this overview, we will concentrate on the application 
of field-theoretic RG techniques to the nonequilibrium dynamics of reaction-diffusion 
systems. Such models consist of classical particles on a lattice, which evolve by hopping 
between sites according to some specified transition probabilities. The particles can also 
interact, either being created or destroyed on a given site following prescribed reaction 
rules (for general reviews, see [SI IE])- However, in order to successfully employ RG 
methods to these systems, one necessary feature is a mean-field theory that is valid 
for some parameter range, typically when the spatial dimension d exceeds an upper 
critical dimension d c . This property renders the renormalization group flows accessible 
within perturbation theory, via a dimensional expansion in e = d c — d. For reaction- 
diffusion models this requirement is usually easily satisfied, since the mean-field theory 
is straightforwardly given by the 'classical' rate equations of chemical kinetics. 

It is the purpose of this topical review to provide an introduction to the methods of 
the field-theoretic RG in reaction-diffusion systems, and to survey the body of work that 
has emerged in this field over the last decade or so. Other theoretical and numerical 
simulation approaches will not be as systematically covered, though various results 
will be mentioned as the context requires. A review by Mattis and Glasser [7j also 
concerns react ion- diffusion models via Doi's representation, but does not address the RG 
methods presented here. Recent overviews by Hinrichsen [HJ and Odor are primarily 
concerned with classifying universality classes in nonequilibrium react ion- diffusion phase 
transitions via Monte Carlo simulations. We will also touch on this topic in our review 
(for brief summaries of the RG approach to this problem, see also Refs. JOJE])- The 
field theory approach to directed and dynamic isotropic percolation, as based on a 
mesoscopic description in terms of Langevin stochastic equations of motion, is discussed 
in depth in Ref. ^2]. As with the previous reviews, our presentation will concentrate 
on theoretical developments. Unfortunately, experiments investigating fluctuations in 
react ion- diffusion systems are deplorably rare. Three notable exceptions are: (i) the 
unambiguous observation of a t -1 / 2 density decay (in an intermediate time window) 
for the diffusion-limited fusion process A + A — > A, as realized in the kinetics of laser- 
induced excitons in quasi one-dimensional N(CH3)4MnCl3 (TMMC) polymer chains |13j ; 
(ii) the demonstration of non-classical A + B — > C kinetics, with an asymptotic t~ 3 / 4 
density decay in three dimensions in a calcium / fluorophore system [Tlj; and (iii) the 
identification of directed percolation critical exponents in studies of spatio-temporal 
intermittency in ferrofluidic spikes [To] . 

The review is organized as follows. The following section provides a basic 
introduction to reaction-diffusion models and the various approaches employed for their 
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investigation. Section |H] describes the mapping of classical reaction-diffusion models onto 
a field theory, while section 0] presents the RG techniques in the context of the kA — > I A 
(£ < k) annihilation reaction. A selection of other single-species reactions is treated in 
section where variations such as Levy flight propagation and the influence of disorder 
are also considered. In addition, this section covers the two-species annihilation reaction 
A + B —>■ with homogeneous and segregated initial conditions, as well as with disorder 
and shear flow. Other multi-species reactions that exhibit similar asymptotic decay are 
also discussed here. Section El deals with directed percolation, branching-annihilating 
random walks, and other examples of nonequilibrium phase transitions between active 
and inactive / absorbing states. The influence of spatial boundaries is also addressed 
here. Finally, in section [?l we give our perspective on open problems for future studies. 

2. Basic Features of Reaction-Diffusion Systems 

2.1. Models 

Our goal is to describe local reactions, of either a creation or annihilation type, for 
which the particles rely on diffusion (or nearest-neighbor hopping) to be brought within 
reaction range. Hence, these processes are often referred to as diffusion-limited reactions. 
Some single-species examples are the pair annihilation reaction A + A — > 0, where '0' 
denotes a chemically inert compound, and coagulation v4+v4 — > A. The diffusive particle 
propagation can be modeled as a continuous- or discrete-time random walk, either on a 
lattice or in the continuum. Reactions occur when particles are within some prescribed 
range; on a lattice, they can also be required to occupy the same lattice site. In such 
systems a single lattice site may be subject to an occupancy restriction (to, say, at 
most n max particles per site) or not, and, of course, the lattice structure can be varied 
(e.g., square or triangular in two dimensions). Computer simulations typically employ 
discrete time random walks, whereas, for example, the analysis of the two-species pair 
annihilation reaction A+B — > by Bramson and Lebowitz (discussed in section l2~4*|) uses 
a continuous time random walk on a lattice with unlimited occupation number, but also 
with an infinite reaction rate so that no site simultaneously contains A and B particles 
[To] . With such a variety of microscopic models to represent a single reaction type, it 
is important to determine which features are universal as opposed to those properties 
that depend on the specific implementation of the processes under consideration. 

The most general single-species reaction-diffusion system can be described by means 
of a set of reaction rules, the ith of which reads kiA — > iiA, with non-negative integers 
ki, £i, and where each process occurs with its own rate or probability per time step. 
Notice that this includes the possibility of reversible reactions, for example 2A «-> A 
as represented by (ki,£i) = (2,1) and (£^,^2) = (1)2). Similarly, directed percolation, 
which defines a broad universality class of nonequilibrium phase transitions between 
active and absorbing states, can be described by the reactions A + A <->• A and A — ► 0, 
where the critical point is reached through tuning the rates of the A — ► (0, 2 A) reactions 
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(see sectional). 

More generic in chemical systems are two-species processes, for example A + B — > 0, 
which requires particles of different types to meet in order for the reaction to occur. The 
different particle species may or may not have the same diffusion constant. A general 
multi-species reaction may be written as J2j kjAj — > J2j^jAj, where Aj labels the jth 
species, and the most general reaction model is then a set of such multi-species processes. 

With such generality available, it is possible to construct both driven and 
relaxational systems. The former case, which includes directed percolation, typically 
comprises both reactions that increase and decrease the particle number. Depending on 
appropriate combinations of the corresponding reaction rates, the ensuing competition 
can, in the thermodynamic limit and at long times, either result in an 'active' state, 
characterized by a finite steady-state density of particles, or a situation that evolves 
towards the empty lattice. For reactions that require at least the presence of a single 
particle, the latter case constitutes an 'inactive' or 'absorbing' state with vanishing 
fluctuations from which the system can never escape. The continuous transition from an 
active to an absorbing state is analogous to a second-order equilibrium phase transition, 
and similarly requires tuning of appropriate reaction rates as control parameters to reach 
the critical region. As in equilibrium, universality of the critical power laws emerges as 
a consequence of the diverging correlation length £, which induces scale invariance and 
independence in the critical regime of microsopic parameters. Alternatively, relaxational 
cases, such as the single-species pair annihilation reaction A + A — > 0, are ultimately 
decaying to an absorbing state: the empty lattice (or a single left-over particle). 
However, here it is the asymptotic decay law that is of interest, and the scaling behavior 
of the correlation functions in the universal regime that is often reached at large values 
of the time variable t. 

2.2. The origin of universality in relaxational reactions 

Reaction-diffusion models provide a rather intuitively accessible explanation for the 
origin of universality. The large-distance properties of random walks are known to 
be universal, depending only on the diffusion constant, a macroscopic quantity. Decay 
processes, such as v4+v4 — > 0, eventually result in the surviving particles being separated 
by large distances, so at late times the probability of a pair of particles diffusing to close 
proximity takes on a universal form. For spatial dimension d < 2 random walks are 
re-entrant, which enables a pair of particles to find each other with probability one, 
even if they are represented by points in a continuum. Therefore, at sufficiently long 
times the effective reaction rate will be governed by the limiting, universal probability of 
the pair diffusing from a large to essentially a zero relative separation. As it turns out, 
this asymptotically universal reaction rate is sufficient to fully determine the entire form 
of the leading density decay power law, that is, both its exponent and the amplitude 
become universal quantities. 

The situation is different in dimensions d > 2: the probability for the pair of 
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particles to come near each other still has a universal form, but even in close proximity 
the reactants actually meet with probability zero if they are described as point particles 
in a continuum. For any reaction to occur, the particles must be given a finite size (or 
equivalently, a finite reaction range), or be put on a lattice. Since the ensuing finite 
effective reaction rate then clearly depends on the existence and on the microscopic 
details of a short- distance (ultraviolet) regulator, universality is (at least partially) lost. 

Consequently, for any two-particle reaction, such as A + A — > (0, A) or A + B — > 0, 
we infer the upper critical dimension to be d c = 2. There is some confusion on this 
point in the literature for the two-species process A + B — > 0, for which sometimes 
d — 4 is claimed to be the upper critical dimension. This is based on the observation 
that for equal initial A and B particle densities the asymptotic density power-law decay 
becomes ~ t~ d / A for d < 4 and t~ l for d > 4. However, this behavior is in fact fully 
exhibited within the framework of the mean-field rate equations, see section 12.31 below. 
Thus d = 4 does not constitute an upper critical dimension in the usual sense (namely, 
that for d < d c fluctuations are crucial and are not adequately captured through mean- 
field approximations). Yet surprisingly, in the two-species pair annihilation process 
there occurs no marked qualitative change at two dimensions for the case of equal 
initial densities. However, the critical dimension d c = 2 strongly impacts the scenario 
with unequal initial densities, where the minority species decays exponentially with a 
presumably nonuniversal rate for all d > 2, exhibits logarithmic corrections in d c = 2, 
and decays according to a stretched exponential law [THl HI] with universal exponent 
and probably also coefficient [TBI EE] for d < 2 (see section fO|) . The critical dimension 
is similarly revealed in the scaling of the reaction zones, which develop when the A and 
B particles are initially segregated (section l5?2|) . 

The upper critical dimension is not always two in diffusion-limited processes. For 
three-particle reactions, such as 3A — > 0, the upper critical dimension becomes d c = 1, 
via the same mechanism as described above: for three point particles to meet in a 
continuous space, they must be constrained to one dimension. For a kth order decay 
reaction kA — > I A (with £ < k), the upper critical dimension is generally found to be 
d c = 2/(k — 1) [20J. Consequently, mean-field descriptions should suffice in any physical 
dimension d > 1 for k > 3. However, this simple argument is not necessarily valid 
once competing particle production processes are present as well. For example, for the 
universality class of directed percolation that describes generic phase transitions from 
active to absorbing states, the upper critical dimension is shifted to d c = 4, as a result of 
combining particle annihilation and branching processes (see section E}. As mentioned 
before, universal features near the transition emerge as a consequence of a diverging 
correlation length, just as for equilibrium critical points. 

Lattice occupation restrictions typically do not affect universality classes for 
relaxational reactions, since the asymptotically low densities essentially satisfy any 
occupation restrictions. A few exceptions are noteworthy: 

(i) The asymptotic decay law in the A + B — > reaction with equal A and B densities 
and in d < 4 dimensions depends on the fluctuations in the initial conditions, which in 
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turn are expected to be sensitive to lattice occupation restrictions [2H [22] • 

(ii) In systems with purely second- or higher-order reactions, site occupation restrictions 
crucially affect the properties of the active phase and the phase transition that separates 
it from the absorbing state [23 EH EH] , see section |H1 

(iii) One-dimensional multi-species systems can be constructed in which the spatial 
ordering of the reaction products specified by the model, which cannot subsequently be 
changed by the dynamics, does affect the asymptotic properties [213 1211 ■ 



2.3. Rate equations 

In general, kinetic rate equations are obtained by taking the rate of change of a given 
species' density or concentration to be proportional to the appropriate product of the 
reactant densities and the reaction rate. This effectively constitutes a factorization of 
higher-order correlation functions (the joint probability of finding a given number of 
reactants at the same location at a given time), and hence corresponds to a mean- 
field type approximation. For example, in the process kA — > £A the probability of a 
reaction is assumed to be proportional to a(t) k , where a{t) denotes the overall (mean) 
A particle density at time t. Such a description that entirely neglects correlations 
and spatial variations is in general justified only if the reactants remain uncorrelated 
and homogeneously distributed (and well-mixed for different participating species) 
throughout the system's temporal evolution. The corresponding rate equation then 
reads 

d t a(t) = -{k-£)Xa{t) k , (1) 

where A represents a reaction rate constant, and the loss rate is proportional to k— £, the 
number of particles removed by the reaction. We assume that £ < k, so d t a is negative. 
Notice that in contrast to k, the integer variable £ does not enter the functional form of 
the rate equation. With an initial density a , Eq. (JTJ is solved by 

a{t) = [l + a k -\k-l)(k-£)\t]^) ' (2) 

which for t ^> l/Xa^ 1 leads to the asymptotic decay a ~ (\t)^ 1 ^ k ~ 1 \ independent of 
the initial density a^. 

Next consider an inhomogeneous system with a local density a(x, t) that is assumed 
to be slowly varying on the scale of the capture radius or lattice size. The rate equation 
approximation for uncorrelated reactants can still be applied; however, the local density 
may now evolve not just via the reactions but also through diffusive particle motion. 
Since the latter process is simply linear in the density, we directly add a diffusion term 
to the rate equation, 

d t a(x, t) = £>V 2 a(x, t) - (k - £)\ a(x, t) k . (3) 

For two-species reactions a pair of rate equations is required. For example, the pair 
annihilation process A + B — > is represented through 

d t a = D A V 2 a - Xab , d t b = D B V 2 b - Xab (4) 
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for the local particle densities a(x, t) and 6(x, t). If both densities may be taken to 
be uniform throughout the temporal evolution, then their decay is just described by 
a(t) ~ b(t) ~ (At) -1 , as in Eq. Q with k = 2. However, this assumption is not 
always justified [21] • Notice that the difference in A and B particle numbers is locally 
conserved by the annihilation reaction. Correspondingly, in the case of equal diffusion 
constants {Da — D b ), the difference field a(x, t) — 6(x, t) satisfies the diffusion equation. 
Thus spatial inhomogeneities relax rather slowly. For the case of equal initial A and B 
densities no, Toussaint and Wilczek (TW) observed that the fluctuations in the initial 
density of this diffusive field in fact determine the long-time characteristics of the two- 
species annihilation process j2H]- With the additional assumption of asymptotic particle 
segregation into separated A and B rich domains, which makes the spatially averaged a 
and b densities exactly half of the \a — b\ average, TW obtained the long-time behavior 

This decay is considerably slower than the one predicted by the uniform rate equation, 
and thus dominates for t — > oo, provided d < 4. Again we emphasize that the 
ensuing qualitative changes in four dimensions merely reflect the importance of spatial 
inhomogeneities within the mean-field rate equation. 



2.4- Relation between RG and other methods 

The rate equations of the previous section still play an important role in the 
field-theoretic analysis that aims to systematically include spatial flucuations and 
correlations. In particular, the rate equation solution represents the zeroth-order term 
in a loop expansion for the density, also called the tree diagram sum. Above the upper 
critical dimension, the higher-order terms in the loop expansion only serve to modify 
the rate constant in some nonuniversal way. Thus it was shown in Ref. [21] that the rate 
equations are asymptotically valid without approximation when d > 2. For d < 2, 
the higher-order terms in a loop expansion provide divergent corrections, which then 
must be regulated (for example, by introducing a lattice) and RG methods brought to 
bear to extract a systematic e expansion. In this case, it is found that the rate equation 
solution gives rise, under RG flow, to the leading-order term in the e expansion. As we 
shall see in subsequent sections, the structure of a loop expansion correcting the rate 
equation solution holds even for more complicated situations, such as the reaction zones 
in the two-species process A + B — > with initially segregated A and B particles, first 
analyzed in Ref. [50] . 

For the A + B — > pair annihilation process, Bramson and Lebowitz demonstrated 
rigorously that the TW decay exponent is exact for all d < 4 for a particular two- 
species model, finding bounds on the density amplitude [H)J H3 EH]. RG methods, 
as mentioned above, showed that the rate equations are asymptotically correct for 
d > 2, and therefore established that the TW result for the density decay, including its 
amplitude, was quantitatively correct and universal for 2 < d < 4. However, while the 
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RG methods suggest that the TW result might extend to d < 2, a demonstration that 
this is indeed the case has not yet, in our opinion, been successfully accomplished (see 
section for details). In this case, it is rather the exact Bramson-Lebowitz result that 
lends credence to the conjecture that the TW decay exponent applies for all dimensions 
d < 4. 

Exact solutions with explicit amplitudes are available for some reaction-diffusion 
models, usually in d — 1. For example, by exploiting a duality with the voter model, 
Bramson and Griffeath solved a particular version of the A + A — > model in one and 
two dimensions i32;, with the result 



The RG calculation for A + A — > generalizes these results for the decay exponent to 
arbitrary dimensions (albeit of limited use for d < 2) [*""""| and provides a quantitative 
expression for the decay amplitude as an expansion in e = d c — d [20] • This expansion 
is poorly convergent in one dimension (e = 1), but, importantly, the RG analysis does 
demonstrate that the density amplitude is universal. Thus, the exact solution and 
the RG analysis contribute in a complementary manner to a full understanding of the 
problem, and tell us that any variation of this model, such as allowing reactions to 
occur whenever reactants are within some fixed number of lattice spacings, will result 
asymptotically in precisely the density decay law Furthermore, the RG expression 
for the amplitude in d = 2 (so e — > 0) is explicit, and matches the exact solution (0). 
This serves both to demonstrate universality of the result and to provide a check on the 
RG method. 

There are many other exact results available in one dimension, and it is beyond the 
scope of this review to provide a complete survey of this field. We restrict ourselves 
to mentioning a few important classes. Many exact techniques exploit a mapping 
from the microscopic master equation onto a quantum spin chain [*""H I*"""*! I37j . 
In several important cases, the ensuing spin system turns out to be integrable and 
a variety of powerful techniques can be brought to bear, such as mapping to free- 
fermion systems jSHl EHl EOJ HTj. These quantum systems may also be studied by 
real space RG methods |42j . The connection between the A + A — > reaction and 
the one-dimensional Glauber dynamic Ising model has also been usefully exploited 
|43 [ 14*4* 1 • F° r steady-state situations, the asymmetric exclusion process, which has 
various reaction-diffusion generalizations, can be solved via a Bethe ansatz or suitable 
matrix ansatz (see Ref. [H"] and references therein). Another useful technique is the 
empty interval method [JZj, which has recently been generalized to a wide range of 
problems |1HJ HHl EH] • Techniques for dimensions d > 1 have also been developed (see 
Ref. [SI] and references therein) . We remark that these mappings generally require that 
each lattice site has finite occupancy (usually zero or one). Although such restrictions 
may originate from physical considerations (e.g., modeling fast on-site reactions) they 
do limit the investigation of universality. Furthermore, very little is currently known 




d= 1 



d = 2 . 



(6) 
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using spin chain mappings about the dynamics of multi-species reactions. 

Besides exact solutions, another important method is Smoluchowski theory [221 EH|, 
which constitutes a type of improved rate equation approximation. Whereas rate 
equations represent a one-point mean-field theory, i.e., a closed equation for the particle 
density, Smoluchowski theory may be viewed as a two-point mean-field approximation, 
namely a closed set of equations for the density as well as the pair correlation function. 
One test particle is taken to be fixed at the origin, and the remaining reactants are 
effectively treated as a non-interacting diffusion field, with the boundary condition 
that their density vanishes at the capture radius of our original particle, and tends 
to some fixed value at infinity. The resulting diffusion flux toward the fixed particle 
is subsequently used to define an effective reaction rate, dependent on the density at 
infinity. The reactive processes are now approximately incorporated by assuming that 
the densities evolve via the usual rate equation but with the new (time-dependent for 
d < 2) effective rate constants. 

This approximation works surprisingly well, in that it actually predicts the correct 
decay exponents for the A + A — > reaction, and even captures the logarithmic 
correction at d c = 2. The amplitudes, however, are incorrect for d < d c , but yield 
reasonable numerical values |5IJ, and are accurate for d = d c J2Hj. Remarkably 
therefore, this improved mean-field theory yields correct scaling exponents even below 
the upper critical dimension. As we shall see in section IH this can be traced to the fact 
that in the corresponding field theory there appears no propagator renormalization, 
and hence no anomalous dimension for the diffusion constant or the fields (see also 
Ref. [SE])- Consequently the density decay exponent turns out to be sufficiently 
constrained (essentially through dimensional analysis) that it is determined exactly, 
i.e., to all orders of the e expansion, within the RG. Yet the density amplitude 
requires an explicit perturbative calculation via the loop (and thus e) expansion. A 
reasonable approximation such as the Smoluchowski theory incorporates the correct 
dimensional analysis that fully determines the decay exponent, but fails quantitatively 
for the amplitude calculation (except at the marginal dimension). Furthermore, mixed 
reactions, such as those considered in section may display decay exponents that are 
not simply fixed by dimensional analysis but rather rely on the details of the particle 
correlations. In these cases, Smoluchowski theory is also insufficient to obtain the 
correct exponents, although here the Smoluchowski exponents have been shown to be 
the same as those from the RG improved tree level |SZj. However, unlike with field- 
theoretic methods, there is no obvious systematic way to improve on Smoluchowski's 
self-consistent approach, with the goal to include higher-order correlations. 
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3. Mapping to Field Theory 

3.1. The model 

We illustrate the mapping to a field theory representation first for the A + A — > single- 
species pair annihilation reaction, and then generalize to other cases. We will consider 
particles on a lattice (say a hypercubic lattice with lattice constant h) performing a 
continuous time random walk, where they hop to a neighboring site at some uniform rate 
D/h 2 (such that D becomes the usual diffusion constant in the continuum limit). The 
particles do not interact, except whenever two or more particles occupy the same site, 
in which case they annihilate with fixed reaction rate A. The state of the system is then 
characterized by the probability P({n}, t) at time t of a particular configuration uniquely 
specified by the string of site occupation numbers {n} = (ni,7i2, •••)■ The system's 
stochastic dynamics is captured through a master equation for the time-dependent 
configuration probability P. For pure diffusion, it assumes the form 

d t P({n}, = E km + -.,^ + 1,^-1...,*) - niP({n}, t) 



h 2 

U <ij> 



+ (nj + 1)P(. . . , ri; - 1, nj + 1 . . . , t) - njP({n}, t) 



(7) 



where the summation extends over pairs of nearest-neighbor sites. The first term in 
the square bracket represents a particle hopping from site i to j, and includes both 
probability flowing into and out of the configuration with site occupation numbers {n} 
as a consequence of the particle move. The second term corresponds to a hop from site 
j to i. The multiplicative factors of n and n + 1 are a result of the particles acting 
independently. 

Combining diffusion with the annihilation reaction gives 

d t P({n},t) = (diffusion term) 
+ \Y^[(n i + 2)(n i + l)P(...,n i + 2,...,t)-n i (n i -l)P({n},t)] , (8) 

i 

where the (diffusion term) denotes the right-hand side of Eq. (JZj). All that remains is to 
specify the initial probability P({n},t = 0). For uniform, random initial conditions the 
particle distribution will be a Poissonian on each site i, i.e., 



no 

Hi 

where n denotes the average number of particles per site. 



nw,o)=nB e - s i (9) 



3.2. Doi's second- quantized representation 

Stochastic classical particle models with local reactions can be rewritten in terms of 
ladder operators familiar from quantum mechanics, as shown by Doi 2J, thus taking 
advantage of the algebraic structure of second quantization. This representation exploits 
the fact that all processes just change the site occupation numbers by an integer. 
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Since we have not implemented any site occupation restrictions, we introduce for each 
lattice site i,j, . . . creation and annihilation operators subject to 'bosonic' commutation 
relations 

[fij, a]] = 5ij , [Si, hj] = [at, aj] = . (10) 

The 'vacuum' (empty lattice) |0) is characterized by Oj|0) = for all i, and on each site 
% we define the state vector \rii) = (a\) ni \0) (note that the normalization differs from the 
standard quantum- mechanical convention). It is then straightforward to show that 

a^ni) = riilrii - 1) , a]\ni) = \rii + 1) . (11) 

Next we employ these on-site vectors to incorporate the state of the entire stochastic 
particle system at time t in the quantity 

i0(*)}=E p (w»*)n(«iri°)- (12) 
{«> » 

As a result, the first-order temporal evolution of the master equation is cast into an 
'imaginary-time Schrodinger equation' 

- d t \<P{t)) = H\<P(t)) , (13) 

where the non-Hermitian time evolution operator ('quasi-Hamiltonian') for the processes 
in Eq. (|HJ) becomes 

H = § £ (3l - - %) - AE[l - («i) 2 ]a? • (14) 

<ij> i 

Equation (|T3J) is formally solved by \4>{t)) = exp(—Ht)\(f)(0)), with the initial state 
determined by Eqs. © and (|T2J). 

The equations of motion for P({n},t) and its moments in this representation are 
of course identical to those following directly from the master equation. Yet at this 
point we may see the advantage of using Doi's formalism: the original master equation 
was complicated by factors of n and n 2 which are now absent. The 'second-quantized' 
formalism provides a natural framework for describing independent particles that may 
be changing in number. 

A different approach is to write down an appropriate Fokker-Planck equation 
[SHI EH HZ] for the processes under consideration, although this is somewhat awkward 
due to the presence of reaction processes where particles are created and / or destroyed. 
However, a more fundamental problem in this approach is encountered in situations 
with low densities. In such cases the Kramers-Moyal expansion used in the derivation 
of the Fokker-Planck equation break down, and it is certainly not valid to terminate 
the expansion in the usual way after the second term [HD]. Alternatively, if the Fokker- 
Planck equation is derived from a coarse-grained Langevin equation, its validity may 
again be questionable, since both the relevant slow variables and their stochastic noise 
correlations must often be inferred phenomenologically. 

We also note that an alternative formalism exists which again starts from a classical 
master equation and leads to a path integral representation. This method uses the 
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Poisson representation, and assumes that the state of the system at time t can be 
expanded into a superposition of multivariate uncorrelated Poissonians [SUl EH 
However, as shown in Ref. [63 , this approach is actually equivalent to Doi's formalism, as 
presented here. An analogous representation in terms of Pauli spin matrices may also be 
used, which then replace the bosonic ladder operators considered here. This corresponds 
to a master equation in which only single occupancy is allowed per particle species at 
each lattice site. These techniques can be especially useful in one dimension, where the 
resulting second-quantized formulation represents certain quantum spin chains, which 
are often integrable j3UE3EEl- However, our primary motive in introducing the second- 
quantized representation here is to map the problem to a field theory, and for this 
purpose the bosonic formalism developed above is more suitable. 

To find ensemble averages A of observables at time t we cannot just use the 
standard quantum- mechanical matrix element, since this would involve two factors of 
the probabilities P. Rather, we need a projection state (V\, defined by the conditions 
(p\a\ = (V\ and (P|0) = 1, leading to 

(p\ = (0|e£i &i . (15) 

From the above properties it follows that 

A(t) = Y,P({n},t)A({n}) =£P({n},t) (V\A\[{^) = {V\A\cj>{t)) , (16) 

where the operator A is given by the function A({n}) with the substitution rij — > ajctj. 
It is naturally not surprising that the expression for averages here differs from usual 
quantum mechanics, since we consider, after all, an entirely classical model. Notice 
too that there is no requirement for H to be Hermitian. In fact, particle annihilation 
and creation reactions obviously lead to non-Hermitian ladder operator combinations. 
Furthermore 

l = (V\<P(t)) = (P\eM-Ht)\m) (17) 

shows that probability conservation is enforced through the condition (V\H = 0. 
Any quasi-Hamiltonian H derived from a probability conserving master equation will 
necessarily satisfy this property. 

The factor exp(J2i a i) in the projection state can be commuted through to the 
right in Eq. (|16p. with the net effect of taking a\ — > 1 + a] [3]. While this has the 
advantage of simplifying the expression for A, we choose not to follow this procedure 
because there are cases, such as branching and annihilating random walks with an even 
number of offspring particles, where such a shift is undesirable as it masks an important 
symmetry. Furthermore, as shown below, the same effect may be obtained in the field 
theory, when desired, by a corresponding field shift. 

We note that probability conservation is reflected in the fact that the time evolution 
operator H must vanish on replacing all a\ — *■ 1. In addition, because of the projection 
state, for any operator A there exists a corresponding operator A', with the same average 
A, that is obtained by normal-ordering A and then replacing the creation operators a\ 
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by unity (the projection state eigenvalue). For example, the density operator d\di may 
be replaced with Oj, and the two-point operator a\di aldj becomes a^j + didj. 

3.3. Coherent state representation and path integrals 

From the second-quantized representation a field theory can be obtained via the very 
same path integral techniques as developed for true quantum many-particle systems. A 
general discussion of these methods can be found in standard textbooks [HUES], and a 
presentation specific to reaction-diffusion models is given by Peliti jlj. For completeness, 
we present the basic method here, with a few supplementary observations. 

First, the (stochastic) temporal evolution is divided into N slices of ultimately 
infinitesimal size At = t/N via 

\<f>(t)) = exp(-tft)|0(O)> = jimexp(-£A£)'/ A *|0(O)> . (18) 

The second-quantized operators are then mapped onto ordinary complex numbers (in 
the bosonic case) or Grassmann variables (for fermions) by inserting a complete set of 
coherent states at each time slice before the limit At — > is taken. 

Coherent states are right eigenstates of the annihilation operator, a|0) = 0|0), 
with complex eigenvalue 0. Explicitly, |0) = exp(— 1|0| 2 + (f>aJ)\0). Their duals are left 
eigenstates of the creation operator: (<fi\a) = (0|0*. A useful overlap relation is 

<0i|0 2 > = exp (-^i| 2 " + <Plfo) ■ (19) 

The coherent states are over-complete, but nonetheless may be used to form a resolution 
of the identity operator. With our convention for the states |n), we have for a single 
lattice site 

i = £ L | n) ( n | = £ L | n) ( m | 5mn =f*± |0) (0| , (20) 

n n - m,n n - J TT 

where we have used 

5 mn = ^- /d 2 0e- |0l V"V\ (21) 
T\m\ J 

with the integration measure d 2 (j) = d(JR(f))d(Q<f)). The above expression generalizes 
straightforwardly to multiple lattice sites according to 

i = |JT^^| W ) (W |, (22) 

where {0} = (0i,02,...) denotes a set of coherent state eigenvalues, one for each 
annihilation operator fij, and |{0}) = |0i) ® |02) <§) — 

A mathematical subtlety can arise with this representation of the identity. 
Expressions such as exp(— HAt) are defined in terms of their power series, with an 
implied sum. The identity operator (|2T)|) contains an integral, and in many cases these 
sums and integrals do not commute. Consider, for example, (0| exp(Aa fc )|0) = 1 with 
the identity (ffilj) inserted: 

£ y^r ! t± mm = £ ^ ffl expH4>m ^ r , (23) 
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where Eq. (J 19)) was used for (O|0). By Eq. (12 lj) the sum over n correctly yields 1+0+0 . . .. 
However, naively exchanging integration and summation gives 

eX p(-|0| 2 + , (24) 



7T 

which does not exist for k > 2 or for k = 2 and |A| > 1. Nevertheless, it is convenient 
to represent expressions such as ()23|) formally by (jUJ), with the understanding that an 
explicit evaluation implies a perturbative expansion in powers of A. As a consequence, 
our formal expressions for the field theory actions will appear to be unstable, but are 
actually well-defined within the framework of perturbation theory. We remark that it 
will be necessary to treat the diffusion terms in H non-perturbatively. This is possible 
because the expansion in powers of the diffusion constant, when acting on the coherent 
state is uniformly convergent in the plane (a sufficient condition for commuting 
sums and integrals), as will be shown below. 

The projection state (V\ is proportional to the dual coherent state with all 
eigenvalues 0j = 1, which for brevity we will call (1|. Also, for Poisson initial conditions, 
|0(O)) is proportional to the coherent state |n ) with all 0j = n - Now label each time 
slice in (II 8j) by a time index r that runs in steps of At from time zero to t and insert 
the identity (j2~2~|) between each slice into Eq. (fTBj) : 

A(t) =A/--Mim J (Y[d%„ T ) (l|i|{0} t ) 

n ({0} T |exp(-#At)|{0} T _ At ) ({0} o |n o ). (25) 

r=At 

The normalization factor Af is to be determined later by requiring the identity operator 
to average to unity. Note that a set of states has been inserted at the time slices at 
and t as well, for reasons that will become clear below. We now proceed to analyze the 
various contributions. 

First, with our caveat above about implied perturbative calculations, we may take 

({0} r |exp(-#At)|{0} r _ A4 ) = ({0} r |{0} T _ A i)exp(-i/({0*} T ,{0} r _ A4 )At) , (26) 

where if({0*} r , {0} T _ A i) = ({0}t|-^1{0}t-a<)- This function is straightforwardly 
obtained by normal ordering H and acting the fij to the right and the d\ to the left, 
whereupon the creation / annihilation operators become respectively replaced with the 
coherent state eigenvalues 0* / 0j. The remaining overlap in Eq. ()26|) factors: 

({<P}r\{<P}r-At) = I] » (27) 

i 

where, according to Eq. |T9|) . for each lattice site 

(4>i t r\4>i,r-At) = eXp(-0* r [0j )T - 4>i, T -At\) exp Q|0i,r| 2 ~ ^\<Pi,r-Atf^j ■ (28) 

Stringing together a product of these states for increasing r will cause the second 
exponential term to cancel except at the initial and final times. The first exponential 
yields a factor exp[— <p ijT *(d(f) ijT /dt)At + 0(At 2 )} for each time slice r and lattice site i. 
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The operator A is assumed to be a function of the annihilation operators a, 
only, by the procedure described at the end of section l3~2l Therefore (l|A|{0} t ) = 
(\\{<j)}t)A({<j)}t), where the latter function is obtained from A through the replacement 
a i — * 4>i,t- The matrix element is multiplied with the remaining exponential factors at 
time t from Eq. (j28|) . giving 

(l|{0}i)nexp Q|0mI 2 ) «exp ($> M ) ■ (29) 

The initial term also picks up a factor from Eq. ()28|) . 

<{0}o|^o> Ilexp ^-^|0 4 ,o| 2 ) oc exp ^ [n o 0* o - |0i,o| 2 ]^) • (30) 

We may now take the limit At — > 0. The O(At) time difference in the 0* and <p T ~At 
arguments of H is dropped with the provision that, in cases where it matters, the 0* 
field should be understood to just follow the field in time. Indeed, this was found to 
be an essential distinction in a numerical calculation of the path integral 66J , and it 
will also play a role in the treatment of the initial conditions below. The O(At) terms 
in the product over r in Eq. (}25|) will become the argument of an exponential: 

A(t)=Af- 1 J (n A({<j>} t ) exp[-S({0*},{0}) o ] • (31) 

i 

Here S denotes the action in the statistical weight, which reads explicitly 

s(m, my = e ( - - mw) + i<mo)i 2 



+ jT' dt[<f>* + #({0*}, {0})] j , (32) 



where we have renamed the final time t — > tf for clarity. The T>cf){D(f)* represent 
functional differentials obtained from Y\ T d(f>i^ T d(f)* T in the limit At — > 0. At last, the 
normalization factor is now fixed via Af = J X\ i T>(\),{D^\ exp[— S({4>*}, {</>})]. 

Before specifying H, let us discuss the initial and final terms in the action ()32j) 
in more detail. The initial terms are of the form exp(0*(O)[0(O) — no}), which implies 
that the functional integral over the variables 0*(O) will create S functions that impose 
the constraints 0(0) = no at each lattic site. Thus the initial terms may be dropped 
from the action (J32|) in lieu of a constraint on the initial values of the fields (f>i [3]. 
However, a path integral with such an implied constraint is not directly amenable to a 
perturbation expansion, so an alternative approach was developed jHZI- All calculations 
will be performed perturbatively with respect to a reference action So composed of 
the bilinear terms oc <$*<$ in S. As will be demonstrated below, any such average will 
give zero unless every factor (p in the quantity to be averaged can be paired up with 
an earlier 0* (that is, the propagator only connects earlier 0* to later 0). The initial 
terms in the action (|32|) can be treated perturbatively by expanding the exponential. 
Recalling that the time ordering of the product 0*(O)0(O) has slightly earlier than 0*, 
we see that all terms in the perturbative expansion will give zero, which is equivalent 
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to simply dropping the 0*0 initial term from the action (|32j) . The remaining initial 
state contribution exp[— n o 0*(O)] then replaces an implied constraint as the means for 
satisfying the assumed random (Poissonian) initial conditions. 

Prior to commenting on the final term — 0;(t/) in the action (J3~2")) . we now proceed 
to take the continuum limit via J2i ~ > I h~ d d d x, <fii(t) — > 0(x, t)h d , and <fi*(t) — » 0(x, t). 
The latter notation indicates that we shall treat the complex conjugate field 0(x, t) and 
0(x, t) as independent variables. This is especially appropriate once we apply a field 
shift 0^1+0, which, in addition to modifying the form of H, has the effect of replacing 

dt dt4> -> <p{t f ) - 0(0) + f tf dt d t <p . (33) 

J 

Thus the final term — 0(x, tf) in the action is cancelled, which simplifies considerably 
the perturbative calculations, but introduces a new initial term. However, the latter will 
again vanish when perturbatively averaged against the bilinear action So, as described 
above. For many of the problems discussed here such a field shift will be employed. 
Lastly, the remaining initial time contribution reads no — > noh d , where no denotes the 
number density per unit volume. Notice that we have (arbitrarily) chosen 0(x, t) to 
have the same scaling dimension as a density. While the continuum limit could have 
been defined differently for the fields and 0, our prescription ensures that the 'bulk' 
contributions to the action must vanish as — > 1 owing to probability conservation. 

At this point, let us explicitly evaluate H for diffusion-limited pair annihilation, 
v4 + v4 — > 0. Since the time evolution operator ()14j) is already normal ordered, we obtain 
directly 

H(m, {0}) = £ £ (fl - m& - fc) - A E(! - 4>f)ti ■ (34) 

11 <ij> i 

We now proceed from a lattice to the continuum limit as outlined above, replacing the 
finite lattice differences in Eq. (jHlj) with spatial gradients. The resulting field theory 
action, prior to any field shift, reads 

S[i, 0] = J d d x |-0(t/) + f * dt [0 (d t - DV 2 ) - A (l - 2 ) 2 ] - n o 0(O)} , (35) 

where Ao = \h d . After applying the field shift — > 1 + 0, we obtain 

S[0, 0] = J d d x {J* f dt [0 (dt - DV 2 ) + A x 00 2 + A 2 2 2 ] - n o 0(O)| , (36) 

with Ai = 2Ao and A 2 = Ao- Finally, we remark again that these actions are defined 
through the perturbation expansion with respect to the nonlinearities, as discussed 
above. However, the diffusion terms are uniformly convergent, resumming to give 
exp(— D\ V0| 2 ), which is bounded for all 0. Hence, the diffusion part of the action 
may be treated non-perturbatively. 



3-4- Generalization to other reactions 



This procedure to represent a classical stochastic master equation in terms of a field 
theory can be straightforwardly generalized to other locally interacting particle systems, 
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e.g., the kth order decay reaction kA — *> £A with £ < k. The appropriate master equation 
for identical particles will result in the time evolution operator 



H = H 



D 



(37) 



where Hp denotes the unaltered diffusion part as in Eq. ()14)) . Following the method 



described above, and performing the field shift 
theory action 



1 + eventually results in the field 



with A,- 



S[$,<j>] = 
Ac (*) - A 




(ft 



- DV 



k 

Ea 

1=1 



-no0(O)j , (38) 
for i > i (note that always 



for i < £, and Aj = A f * 
Afe = Ao). Also, the integer k determines which vertices are present, while £ only modifies 
coefficients. In the simplest case, k = 2, we recover Ai = 2Ao for pair annihilation 
A + A — > 0, whereas Ai = A for pair coagulation A + A — > A. One variant on the 
A + A-^0 reaction would be to allow for mixed pair annihilation and coagulation. That 
is, whenever two A particles meet, with some probability they annihilate according to 
v4 + v4 — > 0, or otherwise coagulate, A + A—>-A. In the master equation these competing 
processes are represented by having both reaction terms present, with reaction rates A^ 
(where £ — 0, 1 indicates the number of reaction products) in the correct proportions. 
The end result is an action of the form (|H8jl above, but with a coupling ratio A1/A2 that 
interpolates between 1 and 2. 

The description of multi-species systems requires, at the level of the master 
equation, additional sets of occupation numbers. For example, the master equation 
for the two-species pair annihilation reaction A + B — > employs a probability 
P({m}, {n},t) where {m}, {n} respectively denote the set of A/B particle occupation 
numbers. Various forms of occupation restrictions could be included in the master 
equation, e.g., Bramson and Lebowitz ^Bj consider a model in which a given site can 
have only A or only B particles. Here we will consider unrestricted site occupation. The 
A and B particles diffuse according to Eq. (jJJ), though possibly with distinct diffusion 
constants. Combining the reactions then gives 



d t P = (diffusion terms) + AE[( m « + + • • > m x + 1j • • • > n i + 1) • • • 1 1) 



(39) 



- miUiP^m}, {n},t) 

The second-quantized formulation then requires distinct creation and annihilation 
operators for each particle species. The state vector is therefore constructed as 

\m= E p({m},{n},t)Y[(dir(bir\o) , m 

{m},{n} i 

and the master equation again assumes the form (|13|). with the time evolution operator 
— ^(a! -a})(a i -a i ) + — 



H 



Effi-^fe-y-AEl 1 

<ij> i 



a\b\)aibi . (41) 
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In the mapping to the field theory we must then involve two sets of coherent states, 
resulting in two independent fields a(x, t) and b(x, t). Hence, after shifting both a — > l+a 
and b — > 1 + b, the action reads: 

(d t -D A V 2 )a + b(d t -D B V 2 )b 



S[a, a, b, b] 




+ A (a + b)ab + A abab 



- a a(0) - 6q6(0) . (42) 



Further generalizations are straightforward: for each new particle species additional 
occupation numbers, second-quantized operators, and fields are to be introduced. The 
details of the reaction are coded into the master equation, though after some practice, 
it is actually easier to directly start with the Doi time evolution operator, as it is a 
more efficient representation. The general result is as follows: For a given reaction, two 
terms appear in the quasi-Hamiltonian (as in the original master equation). The first 
contribution, which is positive, contains both an annihilation and creation operator for 
each reactant, normal-ordered. For example, for the A+A — > and v4+A — > A reactions 
this term reads a^a 2 , whereas one obtains for the A + B — > reaction a^b^ab. These 
contributions indicate that the respective second-order processes contain the particle 
density products a 2 and ab in the corresponding classical rate equations. The second 
term in the quasi-Hamiltonian, which is negative, entails an annihilation operator for 
every reactant and a creation operator for every product, normal-ordered. For example, 
in A + A — > this term would be a 2 , whereas for A + A — > A it becomes a) a 2 , and 
for A + B + C — > A + B it would read crftfabc. These terms thus directly reflect the 
occurring annihilation and creation processes in second-quantized language. 



3. 5. Relation to stochastic partial differential equations 

In some cases, the field theory developed above can be cast into a form reminiscent of 
stochastic partial differential equations (SPDE) with multiplicative noise. Consider the 
action (j3*U|) for single-species pair annihilation A + A — * 0: apart from the quartic term 
A 2 2 2 every term in S is linear in the field. The quartic term can in fact also be 
'linearized' by means of introducing an auxiliary field, where 

exp(— A 2 2 2 ) oc J dr] exp(—r] 2 /2) exp(ir] \J2\ 2 00) . (43) 

Substituting this relation into the action results in three fluctuating fields, namely 0, 
0, and 7], but with the benefit that appears just linearly. Therefore, performing the 
functional integral / T>cf) exp(0[. . .]) simply yields a functional Dirac 5 function, ensuring 
that all configurations 0(x, t) satisfy the corresponding constraint given by its argument. 
As a result, the field is determined by a stochastic partial differential equation: 

d t <j> = D\/ 2 (p-2\ ^ 2 + i^2\' (f)r] , (44) 

where r] represents a stochastic Gaussian variable with unit variance, i.e., (rj) = 0, 
(r)(x.,t) r](x.',t')) = <5(x — x') 5(t — t'). Notice that the above procedure is just the reverse 
of the standard field theory representation of a Langevin-type SPDE pQ. 
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It is not surprising that the noise term vanishes as the field — > 0, since the shot 
noise should diminish as the number of particles decreases. Once the absorbing state 
with zero particles is reached, all deterministic as well as stochastic kinetics ceases. 
However, the appearance of the imaginary noise appears rather strange, and even 
when the stochastic force is redefined as ( = 2^/2Ao0?7, its variance becomes formally 
negative: (C(x, t) C(x', t')) = — 2A 0(x, t) 2 £(x — x') bit — t'). Recall, however, that the 
field is complex, and therefore cannot be simply interpreted as the particle density. 
Furthermore, there is evidence suggesting that the SPDE f)44j) is numerically unstable 

mi eh. 

Interestingly, though, one may obtain an SPDE for a real density field as follows 
|69j . Starting with the unshifted action ([35)1 . we apply the nonlinear Cole-Hopf 
transformation = e p , = pe~?, such that 00 = p, and where the Jacobian is 
unity. This yields 0<9 4 = d t [p(l — p)] + pd t p, and, omitting boundary contributions, 
-D0V 2 = -D0V 2 = -Dp[V 2 p + (Vp) 2 ] for the diffusion term. Finally, the 
annihilation reaction is represented through — Ao(l — 2 )0 2 = Aop 2 (l — e~ 2 ^) = 
2A pp 2 — 2A p 2 p 2 . . ., if we expand the exponential. Hence we see that the quadratic 
term in the field p is now of the opposite sign to before, and therefore corresponds to 
real, rather than imaginary, noise. However, the truncation of this expansion at second 
order is not justifiable, and, furthermore, a consistent description of the annihilation 
kinetics in terms of the fields p and p comes with the price of having to incorporate 
'diffusion noise', i.e. the nonlinear coupling —Dp(Vp) 2 . 

At any rate, this analysis and the previous discussions show that simply writing 
down a mean-field rate equation for the annihilation reaction and then adding real 
Gaussian noise does not in general yield an appropriate SPDE. An even more significant 
observation is that only two-particle reactions can straightforwardly be cast in the form 
of an SPDE, since the linearization (J4"3*j) requires that the field appears quadratically. 
For example, the triplet annihilation reaction 3A — > reaction cannot simply be 
represented as the corresponding rate equation plus real, multiplicative noise. 

4. Renormalization Group Method 

In this section, we describe the basic methodology for performing perturbative RG 
calculations in the context of reaction-diffusion field theories. While our aim is to 
present these techniques in a pedagogical manner, our discussion cannot be entirely 
self-contained here. For additional details, specifically with respect to perturbation 
theory and its representation in terms of Feynman diagrams, reference should also be 
made to the standard field theory literature [101 EH E2] . 

4-1. Diagrammatic expansion 

The diagrammatic expansion for performing field theory calculations is constructed in 
the standard way: the part of the action that is bilinear in the fields is identified as a 
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Gaussian reference action Sq. All other terms are evaluated perturbatively by expanding 
the exponential exp(— S + S ) and averaging with statistical weight exp(— So). These 
Gaussian averages decompose into products of pair correlation functions, and can be 
represented symbolically through Feynman diagrams. Propagators, represented by lines 
in the Feynman graphs, correspond to the pair correlators of fields that are averaged 
together. The nonlinear couplings are graphically depicted as vertices that connect 
propagators together. 

We illustrate this procedure explicitly for the kth order single-species annihilation 
reactions kA — > IA with i < k. In this case, we may take the diffusive part of the 
action as Sq. Thus, the propagator becomes just the diffusion Green function. For a 
single-species reaction, the diffusion action reads 

S = J d d x dt (d t - DV 2 ) = j 0(-p, -co) (-ico + Dp 2 ) 0(p, u) , (45) 

where we have used the time and space Fourier transformed fields 

0(p, lu) = J d d x J dt exp(— ip ■ x + iut) 0(x, t) , (46) 

and extended the time integration range to the entire real axis (as will be justified below) . 
Next we define the propagator as G(x, t) = (0(x, i)0(O, 0)) . Its Fourier transform has 
the form (0(p, u)4>(p', u/))o = Gq(p, oj) {2n) d 5{p + p') 27r<5(a> + u'), with the S functions 
originating from spatial and temporal translation invariance. Explicitly, we infer from 
the action (J4~5|) 

G (p,u)= . ] n 2 , (47) 
— too + Dp z 

which follows from straightforward Gaussian integration. In the complex frequency 
plane the function Gq(p,uj) has a single pole at ou — —iDp 2 . Upon performing the 
inverse temporal Fourier transform from uj to t we find 

G (p,t) =ex P (-Dp 2 t) e(t) , (48) 

where 0(t) denotes Heaviside's step function. Mathematically, its origin is that the sign 
of t determines whether the integration contour is to be closed in the upper or lower 
frequency half plane. Physically, it expresses causality: the unidirectional propagator 
only connects earlier (f> fields to later <p fields (as advertised before in section 13. 3|) . Since 
there exists no earlier source of 4> fields, the time integral in Sq may as well be extended 
from [0, if] to all times, as claimed above. Obviously, for multi-species systems there is 
a distinct propagator of the form (jUj) , (|4l?jl for each particle type. 

The vertices in the Feynman diagrams originate from the perturbative expansion 
of exp(— S + Sq). For example, the term \i<j) l (j) k in the action (|3*K|) requires k mcoming 
propagators, that is, averages connecting k earlier fields to the k later 'legs' attached 
to the vertex, and i outgoing propagators, each of which links a to a later field. The 
diagrammatic representation of the propagator and some vertices is depicted in figure 
Propagators attach to vertices as distinguishable objects, which implies that a given 
diagram will come with a multiplicative combinatorial factor counting the number of 
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time 



-x n 
(a) 




-X-i 



(b) 



-X, 



(c) 




fd) 



Figure 1. Various diagrammatic components: (a) the propagator Go and the initial 
density no, (b) vertices for 2 A — ► £A, (c) vertices for 3 A — > £A, and (d) vertices 
for branching reactions, such as A — > 2A. Our convention throughout is that time 
increases to the left. 



ways to make the attachments. (Note that we do not follow the convention of defining 
appropriate factorials with the nonlinear couplings in the action to partially account 
for this attachment combinatorics.) Any contribution of order m in the coupling Aj is 
represented by a Feynman graph with m corresponding vertices. Moreover, if we are 
interested in the perturbation expansion for cumulants only, we merely need to consider 
fully connected Feynman diagrams, whose vertices are all linked through propagator 
lines. Lastly, the so-called vertex functions are given in terms of one-particle irreducible 
diagrams that do not separate into disjoint subgraphs if one propagator line is 'cut'. 

The perturbative expansion of the initial state contribution exp[— no4>(t = 0)] 
creates <f> fields at t — 0. A term of order n™ comes with a factor of l/m\, but will 
have m\ different ways to connect the initial (f) fields to the corresponding Feynman 
graph, so the end result is that these factorials always cancel for the initial density. A 
similar cancellation of factorials happens for the vertices; for example, a term of order 
A™ does have the 1/m! cancelled by the number of permutations of the m Ai vertices in 
the diagram. Notice, however, these combinatorial factors are distinct from those that 
arise from the different possibilities to attach propagators. 

Since the systems of interest are frequently translationally invariant (in space and 
time), the mathematical expressions represented by the Feynman graphs are often most 
conveniently evaluated in Fourier space. To calculate, say, the mean particle density 
(4>{t)) according to Eq. (pTT]) . one needs diagrams with a single <p field at time t which 
terminates the graph on the left. All diagrams that end in a single propagator line will 
thus contribute to the density, see figure El Since (0(t)) is spatially uniform, the final 
propagator must have p = 0. Similarly, in momentum space the initial density terms 
are of the form n o 0(p = 0, t = 0), so the propagators connected to these also come with 
zero wavevector. The Aj vertices are to be integrated over position space, which creates 
a wavevector-conserving 5 function. Diagrams that contain loops may have 'internal' 
propagators with p ^ 0, but momentum conservation must be satisfied at each vertex. 
These internal wavevectors are then to be integrated over, as are the internal time or 
frequency arguments. 

To illustrate this procedure, consider the second graph in the first row, and the first 
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Figure 2. Feynman graphs that contribute to the mean particle density for the pair 
annihilation and coagulation reactions A + A — > and A + A — > A. The first row 
depicts tree diagrams, the second row one-loop diagrams, and the third row two-loop 
diagrams. 



diagram in the second row of figure El to whose loop we assign the internal momentum 
label p: 

dU G (0, t - h) (-Ai) G (0, hf nl , (49) 



'02 



'12 



t dt 2 [** dt 1 G (0,t-t 2 ) (-Xi, 
Jo Jo 



x J 0-; 2 G (p, h - h) G (-p, h - h) (-A 2 ) G Q (Q, hf n 2 Q (50) 

(the indices here refer to the number of loops and the factors of initial densities 
involved, respectively). The factor 2 in the second contribution originates from the 
number of distinguishable ways to attach the propagators within the loop. Noting that 
Go(p = 0,t > 0) = 1, and that the required p integrations are over Gaussians, these 
expressions are clearly straightforward to evaluate, 

I - XnH I - 8AlA2 "° t2 ~ d/2 (51) 

J ° 2 " ~ Xin ° h2 ~ JMDW 1 (2-d)(4-d) (51) 

(for d 7^ 2,4). Consequently, the effective dimensionless coupling associated with the 
loop in the second diagram is proportional to (\2/D d / 2 )t 1 ~ d / 2 . 

Hence, in low dimensions d < 2, the perturbation expansion is benign at small 
times, but becomes ill-defined as t —>■ oo, whereas the converse is true for d > 2. In 
two dimensions, the effective coupling diverges as (A2/-D) ln(Dt) for both t —>■ and 
t — > 00. The 'ultraviolet' divergences for d > 2 in the short-time regime are easily cured 
by introducing a short-distance cutoff in the wavevector integrals. This is physically 
reasonable since such a cutoff was anyhow originally present in the form of the lattice 
spacing (or particle capture radius). The fluctuation contributions will then explicitly 
depend on this cutoff scale. Thus, in dimensions d > d c — 2, perturbation theory is 
applicable in the asymptotic limit; this implies that the overall scaling behavior of the 
parameters of the theory cannot be affected by the analytic loop corrections, which can 
merely modify amplitudes. In contrast, the physically relevant 'infrared' divergences 
(in the long-time, long-distance limit) in low dimensions d < 2 are more serious and 
render a 'naive' perturbation series meaningless. However, as will be explained in the 
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following subsections, via exploiting scale invariance and the exact structure of the 
renormalization group, one may nevertheless extract fluctuation-corrected power laws 
by means of the perturbation expansion. 




Figure 3. Graphical representation of the Dyson equation for the particle density in 
the A + A — > and A + A — > A reactions. 



Feynman graphs that contain no loops are called tree diagrams. For the particle 
density calculation illustrated in figure 121 these tree diagrams are formed with only Ai 
and n vertices, and we denote the sum of all those tree contributions by a tr (t). For 
example, for the A + A — > (0, A) pair reactions, we may construct this entire series 
iteratively to all orders as shown graphically in figure H3 Thus we arrive at a self- 
consistent Dyson equation for the particle density. More generally, for the single-species 
reactions kA — > IA with i < k the vertex on the right-hand side is connected to k 
full tree density lines. Since all propagators in tree diagrams come with p = 0, the 
corresponding analytical expression for the Dyson equation reads 

a tr (t) = n - \, f dt' a tr (t') k . (52) 
Jo 

Upon taking a time derivative, this reduces to the mean-field rate equation with 
the correct initial condition. Furthermore, Ai = (k — ^)Ao, i.e., the rate constant is 
properly proportional to the number of particles removed by the reaction. Evidently, 
therefore, the tree-level approximation is equivalent to simple mean-field theory, and 
any fluctuation corrections to the rate equation must emerge from Feynman graphs 
that incorporate higher-order vertices Aj with i > 1, i.e., diagrams with loops. We note 
that the mean- field rate equations also follow from the stationarity conditions, i.e., the 
'classical field equations', for the action S (regardless of performing any field shifts). For 
example, taking 5S/S(p = = 5S/5(j) for the action (jHHjl results in = and (ft = a tT (t). 

Following up on our earlier discussion, we realize that the loop fluctuation 
contributions cannot alter the asymptotic power laws that follow from the tree diagrams 
in sufficiently large dimensions d > d c , where mean-field theory should therefore yield 
accurate scaling exponents. However, note that, for d > d c , there will generally 
be non-negligible and non- universal (depending on the ultraviolet cutoff) fluctuation 
corrections to the amplitudes. Recall that the (upper) critical dimension d c can be 
readily determined as the dimension where the effective coupling associated with loop 
integrals becomes dimensionless. 

4-2. Renormalization 

As we have seen in the above example (|5()|) . when one naively tries to extend the 
diagrammatic expansion beyond the tree contributions to include the corrections due 
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to loop diagrams, one encounters divergent integrals. We are specifically interested in 
the situation at low dimensions d < d c : here the infrared (IR) singularities (apparent as 
divergences as external wavevectors p — > and either t — > oo or uo — > 0) emerging in the 
loop expansion indicate substantial deviations from the mean-field predictions. Our goal 
is to extract the correct asymptotic power laws associated with these 'physical' infrared 
singularities in the particle density and other correlation functions. To this end, we shall 
turn to our advantage the fact that power laws reflect an underlying scale invariance 
in the system. Once we have found a reliable method to determine the behavior of any 
correlation function under either length, momentum, or time scale transformations, we 
can readily exploit this to construct appropriate scaling laws. 

There exist well-developed tools for the investigation and subsequent renormaliza- 
tion of ultraviolet (UV) singularities, which stem from the large wavenumber contribu- 
tion to the loop integrals. In our models, these divergences are superficial, since we can 
always reinstate short-distance cutoffs corresponding to microscopic lattice spacings. 
However, any such regularization procedure introduces an explicit dependence on the 
associated regularization scale. Since it does not employ any UV cutoff, dimensional 
regularization is especially useful in higher-loop calculations. Yet even then, in order 
to avoid the IR singularities, one must evaluate the integrals at some finite momentum, 
frequency, or time scale. In the following, we shall denote this normalization momentum 
scale as k, associated with a length scale k~ x , or, assuming purely diffusive propaga- 
tion, time scale t = 1/(Dk 2 ). Once the theory has been rendered finite with respect 
to the UV singularities via the renormalization procedure, we can subsequently extract 
the dependence of the relevant renormalized model parameters on k. This is formally 
achieved by means of the Callan-Symanzik RG flow equations. Precisely in a regime 
where scale invariance holds, i.e., in the vicinity of a RG fixed point, the ensuing ultra- 
violet scaling properties also yield the desired algebraic behavior in the infrared. (For 
a more elaborate discussion of the connections between UV and IR singularities, see 
Ref. H2|.) 

The renormalization procedure itself is, in essence, a resummation of the naive, 
strongly cutoff-dependent loop expansion that is subsequently well-behaved as the 
ultraviolet regulator is removed. Technically, one defines renormalized effective 
parameters in the theory that formally absorb the ultraviolet poles. When such a 
procedure is possible — i.e., when the field theory is 'renormalizable', which means 
only a finite number of renormalized parameters need to be introduced - one 
obtains in this way a unique continuum limit. Examining the RG flow of the scale- 
dependent parameters of the renormalized theory, one encounters universality in the 
vicinity of an IR-stable fixed point: there the theory on large length and time scales 
becomes independent of microscopic details. The preceding procedure is usually only 
quantitatively tractable at the lowest dimension that gives UV singularities, i.e., the 
upper critical dimension d c , which is also the highest dimension where IR divergences 
appear. In order to obtain the infrared scaling behavior in lower dimensions d < d c , we 
must at least initially resort to a perturbational treatment with respect to the marginal 
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couplings in the theory, which are dimensionless at d c (and perhaps subsequently resum 
the perturbation series). The scaling exponents can then be obtained in a controlled 
manner in a dimensional expansion with respect to the small parameter e = d c — d. 

We can follow standard procedures (see, e.g., Refs. [7UJ [7IJ |?2|) for implementing 
the renormalization program. First, we must identify the primitive UV divergences, 
the sub-components of the diagrams that are responsible for these singularities. This 
is most conveniently done using the vertex functions T^ m ' n \ which represent a sum of 
all possible one-particle irreducible Feynman graphs that are attached to n incoming 
and m outgoing propagator lines (of course, for a multispecies vertex function, separate 
indices are required for the incoming and outgoing lines of each species). In frequency 
and wavevector space these subdiagrams enter multiplicatively, which means that once 
the vertex function divergences are resolved, the general diagrammatic expansion will 
be well-behaved. Which vertex functions are primitively divergent can be ascertained 
by direct power counting: [r (m ' n )] = n a , where k denotes some reference wavevector, 
such as the normalization scale mentioned above. The scaling dimension of the vertex 
function T <ym,T ^ is just that of the coupling X mn from a term A mn m n in the action, 
since at the tree level T^ m ' n ^ ~ A mn . Loop diagram fluctuation corrections, however, 
require additional nonlinear couplings, which in turn determine the primitive degree of 
divergence for the associated momentum space integrals. We refer to the standard field 
theory texts for a general discussion of the ensuing vertex function dimensional analysis, 
but provide a brief outline of the procedure for our situation. 

The action S itself must be dimensionless, so any term in the integrand of S 
must have scaling dimension K d+2 . Consider the kA — > P. A annihilation reaction 
with i < k. The interaction vertices Aj</>*0 fc in the action (jHEj) correspond to k 
incoming and i = 1 . . . k outgoing lines. With our choice of taking the continuum 
limit, the scaling dimension of the fields are [(f)] = n d and [(f)] = n°. Hence we obtain 
[p(i,fc)] _ ^.j _ K 2~{k~i)d f Qr a rj ^ (recall that Aj oc A of the unshifted theory). Let 
us now investigate the lowest-order loop correction to the vertex function T^' k \ which 
contains k internal propagator lines, and thus is proportional to AjAfe. The involved 
momentum integral therefore must scale as [AJ -1 = k^" 1 ^ -2 . By choosing as k the 
inverse short-distance cutoff, we see that if the exponent here is non-negative, the vertex 
function contains a primitive UV divergence (as k — > oo) and must be renormalized. 
(The converse is true for the IR singularities, which emerge in the limit k — > 0.) The /-th 
order corrections to the bare vertex function must scale as K l i( k - 1 ) d ~ 2 } . Consequently, for 
a given k, the vertex functions T^'^ become primitively UV-divergent, and increasingly 
so in higher loop orders, for d > d c = 2/(k — 1) for all i. The IR singularities, on the 
other hand, become successively worse in higher orders of the perturbation expansion 
for d < d c . At the critical dimension, the loop diagrams carry logarithmic UV and IR 
divergences, independent of the loop order. 

Hence this scaling discussion already reveals the upper critical dimension d c = 2 
for pair annihilation and coagulation, A + A — > (0, A), in agreement with our analysis 
following Eq. (JoTj) . whereas d c = 1 for the triplet reactions 3 A — > [A. All higher-order 
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reactions should be adequately described by mean-field theory, as represented by the 
tree Feynman graphs. But in general, since the assignment of scaling dimensions to the 
fields and is somewhat arbitrary, one needs to be more careful and first determine 
the effective couplings in the perturbational expansion and from there infer the upper 
critical dimension. For example, the field theory for directed percolation, see section |H1 
below, incorporates the vertices \i 2 4><f) 2 and \ 2 i<j) 2 <j). The effective coupling then turns 
out to be the product u ~ Ai 2 A 2 i, whose scaling dimension is [it] = K 4 ~ d , which indicates 
that actually d c = 4 in this case. 

Once the primitive divergences are identified, they are used to define renormalized 
coupling constants into which the UV singularities are absorbed. For the kA — > I A decay 
reaction, this procedure is unusually simple, and in fact the entire perturbation series 
can be summed to all orders. First, we note that since all vertices in the action (}38|) have 
k > 2 incoming lines, one cannot construct any loop diagram for the vertex function 
r^ 1,1 ) that corresponds to the inverse propagator. Hence the bare diffusion propagator 
(JUj) or (|4*K|) is not affected by fluctuations and the tree-level scaling x ~ p^ 1 ~ (Dt) 1 / 2 
remains intact [33]. The absence of field and diffusion constant renormalization in 
this case constitutes the fundamental reason that improved mean-field theories of the 
Smoluchowski type, and also simple scaling approaches, are capable of obtaining correct 
density decay exponents. Hence, for pair annihilation or coagulation, one would predict 
that the density of surviving particles at time t is given in terms of the diffusion length 
by a(t) ~ x- d ~ (Dt)- d / 2 . 

t 2 t x ^ CXZX 

Figure 4. The set of primitively divergent diagrams contributing to I^ 1 ' 3 )^ — ii) for 
the 3A -> IA reaction, I < 2. 

This leaves us with the renormalization of the vertex couplings A m , encoded in the 
vertex functions r( m ' k \ Since all A m are proportional to the reaction rate Ao, there is 
essentially only a single renormalized coupling gn (as is obvious if we work with the 
unshifted action). As a representative example, we depict the ensuing diagrammatic 
expansion for the case m = 1 and k = 3 in figure 0] It will be sufficient to work with 
p = 0. In (p, t) space we obtain explicitly 

Y^ k \t 2 - t x ) = X m S(t 2 - h) - X m X I(t 2 - h) 

+ A m A 2 f 2 dt' I(t 2 - t')I{t' - ti) - . . . , (53) 

where we have used A& = Ao, and I(t) is given by the following integral with k — 1 loops 
over k propagators: 

= kljfl (j^ {2*) d 8 (X>) exp(-E^) = B k (Dt)-^ , (54) 
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with B k = k\k d / 2 (47r) d l dc and d c — 2/{k— 1). Taking the Laplace transform 

~ poo 

r (ra ' fe) (s) = T^ k \t)e- st dt , (55) 
Jo 

the convolution theorem renders (|53|) into a geometric sum: 

P K*)( s ) = *™ , = 5fc r(e/d c )D- d /^ S -^ , (56) 

1 + X Q I{s) 

where e = d c — d, and T(x) is Euler's gamma function. For e > 0, the loop integral I(s) 
displays the expected IR divergence as s — > 0. On the other hand, the UV singularity 
(for finite s) at d c appears as an e pole in the gamma function T(e/d c ) oc d c /e. 

Recall that the bare effective coupling Xq/D has scaling dimension K 2 -( fc - 1 ) d = 
K 2e l dc . Hence we may introduce the dimensionless parameter g = (Xq/D) n~ 2t l dc with 
some arbitrary momentum scale k. Next we define its renormalized counterpart by 
the value of the vertex function r( fc ' fc ) at vanishing external wavevectors and Laplace 
argument s = 1/to = Dk 2 , which sets our normalization point. Thus, 

9R = (s) \ s=Dk2 k-^/D = Z g9o , Z g l = l + g B k T{e/d c ) (57) 

according to Eq. (JSTIj) . Formally, a perturbative expansion in terms of go can be readily 
exchanged for an expansion in qr by straight substitution g = gn/[l — gnBkT(e/d c )}. 
As e — > 0, however, this substitution becomes singular, since the multiplicative 
renormalization constant Z~ l that was introduced to absorb the UV pole diverges in 
this limit. It is crucial to note that the renormalized coupling explicitly depends on 
the normalization scale k. Indeed, this is borne out by calculating the associated RG (3 
function 

d_ 

dn 

Notice that /3 g , here computed to all orders in perturbation theory, is regular as e — » 
when expressed in terms of renormalized quantities. For the simple annihilation models, 
the P function is exactly quadratic in g R , which is not typical and is due to the 
geometric sum in the primitively divergent vertex function, which reduces the fluctuation 
contributions effectively to the one- loop graph. Eq. (|55|) has the standard structure that 
the linear coefficient is of order e = d c — d while the quadratic coefficient is order unity. 
The theory is manifestly scale-invariant (independent of k) at either gn = (which here 
corresponds to pure diffusion, no reactions) or at the special fixed-point value for the 
coupling given by the nontrivial zero of the j3 function 

g R = [B k T(e/d c )]- 1 , (59) 

which is of order e. 

We remark that the above renormalization structure directly applies to multispecies 
annihilation reactions as well. Consider, for example, the action (pl2*|) for A + B — > 0. 
Once again, the vertices permit no propagator renormalization, and the perturbation 
expansion for the vertex functions that define the renormalized reaction rate takes the 
same form as in figure EJ with incoming A and B lines, and the series of internal loops 
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formed with precisely these two distinct propagators. Consequently, we just recover the 
previous results (|57jl. (JSHJ, and (jo^j) with k = 2. 



^.5. Callan-Symanzik equation and loop expansion for the density 

In order to employ the renormalization group machinery to obtain asymptotic (long- 
time, large-distance) expressions for the particle density and its correlations, we next 
develop the Callan-Symanzik equation for the density. By means of perturbation theory 
in the IR-fmite regime and the subsequent substitutions of g R for A , the density can 
be calculated as an explicit function of time t, initial density no, diffusion constant D, 
renormalized coupling g R , and arbitrary normalization point k (or the equivalent time 
scale to — 1/-Dk 2 ). Since k does not appear at all in the unrenormalized theory, we 
must have n da{t, no, D, Ao) jdn = for the particle density at fixed bare parameters, or 
equivalently, after rewriting in terms of the renormalized density, 

d d 

a(t,n ,D,K,g R ) = (60) 



with the (3 function Yet dimensional analysis tells us that a(t, no, D, k, g R ) = 

K d a(t/to, no/n d , g R ), reducing the number of independent variables to three. 



Consequently, Eq. (|6T)j) yields the Callan-Symanzik (CS) equation 

d ,d „ , , d 



2Dt 



dn t: h P g {g R ) 



+ d 



a(t,n ,D, K,g R ) = 



(61) 



d(Dt) "dn dg R 

This partial differential equation is solved by the standard method of characteristics, 
where one introduces a flow parameter via k — > k£. The IR asymptotic region then 
corresponds to i — > 0. For our purposes, it is most convenient to directly employ 
( K £) 2 = l/Dt or £ 2 = t /t. Thereby we find 

a(t, n , t , g R ) = (t /t) d/2 a(n (t),g R (t)) , (62) 
with the running initial density 

n Q (t) = (t/to) d/2 n , (63) 
and the running coupling g R defined by the solution of the characteristic equation 



d~g R {t) 



-2t 



dg R {t) 



R 



gR{Z = i) = g R (t = t 



9r 



(64) 



dt dt 

The method of characteristics requires a known value of the function, in this case the 
density, for some value of the running parameters. Since we chose the normalization 
point s = l/t > outside the IR-singular region, we may use perturbation theory to 
calculate the right-hand side of Eq. (jo^j) . The Callan-Symanzik equation then allows us 
to transport this result into the perturbatively inaccessible asymptotic region. 

Because of the simple form of the (3 function (p)%|) . the running coupling can be 
found exactly by integrating the flow equation (JHlj): 



9R{t) = 9*r 



1 + 



9r~9r 
9r 



t/d c 



-1 



(6^0). 



(65) 
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Thus we see for e > that the running renormalized coupling approaches the RG fixed 
point (joTJj) as t — > oo, independent of its initial value g R . Thus, universal behavior 
emerges in the asymptotic regime, whose scaling properties are governed by the IR- 
stable fixed point g* R . Therefore an expansion in powers of go is converted, via the CS 
equation, to an expansion in powers of e. For this purpose, we merely need to invert 
Eq. (|SZJ to find g in terms of g R : g = g R /[l - g R /g* R ) = g R + g 2 R /g* R + . . .. Notice that 
9R = 9r formally corresponds to g ~ Xq/D = oo, i.e., the annihilation reactions are 
indeed diffusion-limited. Above the critical dimension (e < 0), g R (t) — > algebraically 
~ t~\ e \l dc 1 whereas precisely at d c the running coupling tends to zero only logarithmically, 

g R (t) = 9R - - - (e = 0) . (66) 

yRK) l + B k g R \n(t/t ) K J 1 ; 

We may use these findings to already make contact with both the rate equation 

and Smoluchowski approximations. For d > d c , the effective reaction rate X(t) ~ 

D(Ki) 2t / dc g R (£) = D(Dt)~ t l dc g R (t) — > const, asymptotically, as implicitly taken for 

granted in mean- field theory. Below the critical dimension, however, X(t ^> t ) ~ 

D{Dt)~ t l dc g R or its density- dependent counterpart X(a) ~ Da 2e ^ dd ^ decrease precisely 

as in the Smoluchowski approach. At d c , we have instead X(t) ~ D/ \n(t/to) or 

X(a) ~ D/ln(l/a). Replacing A — > X(t) or X(a) in the mean-field rate equations (|T|1 

then immediately yields the results ( fB^l for k — 2, whereas a{t) ~ \hi(Dt) / Dt} 1 ' 2 for 

k — 3 at d c — 1. 



(a) (b) 



(c) 



"<o< ^ ^ 

(d) (e) (f) (g) 

Figure 5. One-loop and two-loop Feynman diagrams for the particle density, shown 
for the pair annihilation reaction k = 2. 



While we have now established a systematic expansion in terms of g R , a perturbative 
calculation in powers of no is not useful, since h (t) diverges for t ^ t , Eq. (jfiHjl . It 
is thus imperative to calculate to all orders in the initial density no. To this end, we 
proceed to group the Feynman graphs for the particle density according to the number 
of closed loops involved. First, we obtain the tree diagrams represented by the Dyson 
equation in figure El When substituted into the right-hand side of Eq. (|6*2*|) the limit 
no — > oo will give a finite result, with leading corrections ~ 1/wo ~ t~ d ^ 2 that vanish 
asymptotically. Explicitly, replacing the bare with flowing renormalized quantities in 
Eqs. (JH2J) and ^ at the RG fixed point gives 

a(t) = TO — ► A M (Dty d ' 2 , (67) 

[1 + n /dc {k -l)(k- £)g* R {Dty/**]**/* 

with universal amplitude Am = [(k — !)(& — (^g*^ 1 ^^^ ■ 
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Next, we consider one- (a) and two-loop (b)-(g) diagrams for the density, depicted 
in figure El (all for the case k = 2, but the generalization to arbitrary k is obvious). In 
order to sum over all powers of no, the propagators in these diagrams are replaced with 
response functions, which include sums over all tree-level dressings. These are depicted 
in figure El along with the Dyson equation they satisfy. For d < d c each vertex coupling 
asymptotically flows to the 0(e) fixed point (J59j) . so the loop expansion corresponds to 
an ordering in successive powers of e = d c — d. However, each order of the expansion, 
under RG flow, comes in with the same t~ d l 2 time dependence. Thus, the loop expansion 
confirms that the exponent is given explicitly by the tree-level result, and provides an 
epsilon expansion for the amplitude of the density decay. 

— = + + 




Figure 6. Response function for k=2 



As mentioned above, the very same renormalizations hold for multi-species 
annihilation reactions, for example the pair process A + B — ► 0. Consequently, in the 
case of unequal initial A and B densities (with ao < bo, say) in dimensions d > d c = 2, 
the mean-field result that the minority species vanishes exponentially a(t) ~ exp(— Xt) is 
recovered, whereas for d < 2 the direct replacement A — > A(i) ~ D(Dt)~ 1+d l 2 correctly 
yields a stretched exponential decay a(t) ~ exp[— const. (Dt) d l 2 }, while at d c — 2 
the process is slowed down only logarithmically, a(t) ~ exp[— const. Dtj ln(Dt)]. The 
asymptotic B particle saturation density is approached with the same time dependences. 
The amplitudes in the exponentials were computed exactly by other means in dimensions 
d < 2 by Blythe and Bray 02]. 

5. Further Applications 

Now that we have established the basic field-theoretic RG machinery necessary to 
systematically compute exponents and amplitudes, we can summarize some more 
sophisticated applications. We deal first with systems without phase transitions, before 
moving on in sectionElto describe reaction-diffusion systems that display nonequilibrium 
phase transitions between active and absorbing states. Our aim in this section will be 
to give a brief outline of the results available using RG methods, rather than to delve 
too deeply into calculational details. 

5.1. Single- species reactions 

• The kA — > C.A reaction with £ < k: 

The RG treatment for the general single-species annihilation reactions kA — >• £A (I < k) 
was explicitly covered in the previous sections. The upper critical dimension of these 
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reactions is d c — 2/(k — 1), and we note in particular that the reactions A + A — > and 
A + A — > A are in the same universality class |33j . We also emphasize again that, for 
d < d c , the amplitudes and exponents are universal, independent of the initial conditions 
(apart from highly specialized initial conditions, such as those in Ref. [73*] . where the 
particles were initially positioned in pairs). 

• A + A — > (0, A) with particle input — ► A: 

Droz and Sasvari [JJ] studied the steady state of the combined A + A — > (0, A) and 
— > A reactions, focusing particularly on how the density scales with J, the particle 
input rate. This process appears as an interaction </</> in the action. Power counting 
gives [J] = n 2+d and straightforward arguments show that for d < d c — 2, and for 
sufficiently small values of J, the density scales as j d /( d + 2 )^ and that the characteristic 
relaxation time behaves as r ~ j~ 2 /( d + 2 ). Finally, these findings were combined to 
reproduce the standard density scaling as t~ d l 2 . Rey and Droz extended this approach 
to provide explicit perturbative calculations of the density scaling function [7""] . 

• Disordered systems: 

Another important variation on these simple reaction-diffusion models is to include 
quenched disorder in the transport. Various models of quenched random velocity fields 
in the A + A — > reaction have been investigated using RG techniques, including 
uncorrelated (Sinai) disorder [Jjj and also long-ranged correlated potential disorder 
[77] , We consider first the case of (weak) Sinai disorder, with velocity correlator 
(v a (x)v l3 (y)) = A5 atl 35(x—y) analyzed by Richardson and Cardy [76 . An effective action 
is found by averaging over this disorder, after which one must renormalize the disorder 
strength and diffusion constant in addition to the reaction rate. Unlike the case of pure 
diffusive transport, it turns out that the amplitude for the asymptotic density decay rate 
as a function of time is nonuniversal for d < 2: n ~ Cdt~ d ^ z , with z = 2 + 2e 2 + 0(e 3 ) 
and e = 2 — d, but where Cd must be nonuniversal on dimensional grounds. It is only 
when rewriting the density as a function of the disorder-averaged diffusion length that a 
universal scaling relation emerges: n ~ Bd(r 2 )~ d ^ 2 , where is universal. Results were 
also obtained for d = 2, where the effects of the uncorrelated disorder are not strong 
(only the amplitude, but not the exponent, of the asymptotic density decay is altered). 
Interestingly, for weak disorder, it was found that the amplitude of the density decay 
is reduced, implying that the effective reaction rate is faster than for the case of purely 
diffusing reactants. Physically, this results from the disorder 'pushing' particles into 
the same region of space, thus speeding up the kinetics. Theoretically, this originates 
in a disorder- induced renormalization of the reaction rate. However, as the disorder is 
increased it was also shown that the reaction rate would then begin to decrease. This 
stems from a disorder-induced renormalization of the diffusion, which works to slow the 
kinetics, i.e. operates in the opposite direction to the disorder- induced reaction rate 
renormalization. 

The related case of long-ranged potential disorder, where the random velocity field 
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can be considered as the gradient of a random potential, was analyzed using RG methods 
in two dimensions by Park and Deem [77]. In this case, rather more drastic effects 
were found, with an altered decay exponent from the case of purely diffusing reactants. 
Physically, this results from the different nature of the disordered landscape [7S] , where 
for long-ranged potential disorder, but not of the Sinai type, deep trapping wells exist 
where, in order to escape the trap, a particle must move in an unfavorable direction. 
Park and Deem employed replicas to analyze the effect of long-ranged disorder, where the 
correlation function of the quenched random potential behaves as 7/A; 2 . They obtained 
that the asymptotic density decay was modified to where 5 is defined in the absence 
of reaction by the anomalous diffusion relation (r(t) 2 ) ~ t 1 ^ 5 . Here, 5 was found to be 
a nonuniversal exponent depending on the strength of the disorder. The amplitude of 
the decay also turned out to be a nonuniversal quantity. 

We also mention that 'superfast' reactivity has been found in d = 2 for the 
A + A —> reaction in a model of turbulent flow with potential disorder [7JJ]. This 
case was also investigated numerically jHH|- RG methods indicated that this regime 
persists in a more realistic time- dependent model for the random velocity field [80J. 
The case of A + A — > also in a time-dependent random velocity field, but generated 
now by a stochastically forced Navier-Stokes equation, was considered in Ref. [ST] . 

• Levy flights in reactive systems: 

Replacing diffusive propagation with long-ranged Levy flights constitutes another 
important modification to the dynamics of reactive systems. Such Levy flights are 
characterized by a probability for a particle's jump length I decaying for large i as 
p ~ Cr d -° . For a < 2 this results in a mean-square displacement in one dimension 
growing as t 1//fJ , faster than the t 1 ^ 2 law of diffusion. Naturally, one expects that altering 
the dynamics of the system in this way will modify the kinetics, as one is effectively 
making the system better mixed with decreasing a. The propagator for Levy processes 
becomes Go(p ; w ) — (—iu+D L p a )~ 1 , meaning that time scales acquire scaling dimension 
hT u rather than k~ 2 (for a < 2). Consequently, power counting for the A + A —>■ (0, A) 
reaction gives [A] = K ff " d , implying that d c = o. Once again only the reaction rate is 
renormalized, which then flows to an 0(e = a—d) fixed point under the RG. Dimensional 
analysis subseqently fixes the asymptotic density decay rate as t~ d ^ a , ioi d < o [82J. 

Note that the upper critical dimension is now a function of the Levy index a. This 
feature has been exploited by Vernon jHS] to compute the density amplitude for the 
A + A — > reaction with Levy flights to first order in e = o — d. a was then set to 
be slightly larger than unity and the behavior of the system was studied numerically 
in d = 1. This ensures that e = o — d is a genuinely small expansion parameter (i.e., 
e <C 1) in the physical dimension d = 1. This contrasts with the case of A + A — > with 
standard diffusion where, in order to access d — 1, e — 2 — d must be set to unity. As we 
have seen, in that situation the expansion for the density amplitude agrees only rather 
poorly with numerics |20j. However, for the Levy flight case, Vernon demonstrated 
that the accuracy of the expansion indeed improves with decreasing e (i.e. decreasing 
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cr towards unity). This ability to vary the value of d c has also been used to probe the 
behavior of directed percolation and branching-annihilating random walks [H^IEUIES], 
see section E] below. We also mention that the reaction A + A — > with Levy flights and 
quenched disorder was studied using RG methods in Ref. jHSl- Finally, the authors of 
Ref. |HZ| used RG techniques to investigate the case of short-ranged diffusion, but with 
long-ranged reactive interactions. 

5.2. Two-species reactions 

• The homogeneous A + B — > reaction: 

The two-species decay reaction is perhaps the most relevant to chemical systems. It is 
also considerably more complicated to analyze, since the A + B —>■ pair annihilation 
process leaves the local density difference field a — b unchanged. This conservation 
law provides a slow mode in the dynamics that is crucial in determining the long-time 
behavior of the system. We consider first the case where the A and B particles are 
initially mixed together throughout the system. If their initial densities are unequal, 
say b > a , the asymptotic dynamics will approach a steady concentration of bo — a of 
B particles, with very few isolated A particles suriving. In this situation, exact results 
indicate an exponentially decaying A particle density for d > 2, logarithmic corrections 
to an exponential in d = 2, and a stretched exponential exp(— c^/t) form for d — 1 
[HH HH1 HH], where c is a constant. As briefly discussed in section 14.31 above, these 
results correspond in the RG framework to the standard renormalization of the reaction 
rate. 

In contrast, when starting from equal initial densities, the fluctuations in the initial 
conditions for the difference field a — b decay to zero slowly, by diffusion. This case was 
studied by Toussaint and Wilczek [2H| based on the idea that after a time t, on length 
scales shorter than the diffusion length Id ~ t 1 ^ 2 only whichever of the species happened 
to be in the majority in that region initially will remain. In other words, the two species 
asymptotically segregate. Since the initial difference between the A and B particle 
numbers in that region is proportional to Z^ 2 , this leads to an asymptotic t _d//4 decay 
|29j . Clearly, for d < 4, this dominates the faster t _1 mean-field density decay which 
assumes well-mixed reactants throughout the system's temporal evolution. Toussaint 
and Wilczek explicitly calculated the amplitude for this decay under the assumption 
that the only relevant fluctuations are those in the initial conditions. These results were 
since confirmed by exact methods [TBJ El EH] • 

Turning to the field-theoretic RG approach, the action (|42|) for the process A + B — > 
contains diffusive propagators for both A and B species, possibly with unequal diffusion 
constants, together with the interaction vertices Xaab, Xbab and Xabab. Power counting 
reveals [A] = K 2 ~ d , the same as in the A + A — > reaction. This implies that the upper 
critical dimension is d c = 2, consistent with the behavior for unequal initial densities. 
The renormalization of the A + B — ► action also follows similarly to the A + A —>■ 
case. Surprisingly, however, a full RG calculation of the asymptotic density in the case 
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of unequal initial densities has not yet been fully carried through. For the equal density 
case, though, a field theory approach by Lee and Cardy is available [21] . The Toussaint- 
Wilczek analysis reveals a qualitative change in the system's behavior in four dimensions, 
whereas the field theory yields d c = 2. The resolution of this issue lies in the derivation 
of an effective theory valid for 2 < d < 4, where one must allow for the generation 
of effective initial (t = 0) 'surface' terms, incorporating the fluctuations of the initial 
state. Aside from this initial fluctuation term, it was shown that the mean-field rate 
equations suffice [21]. Using the field theory approach, Lee and Cardy were also able 
to demonstrate the asymptotic segregation of the A and B species, and thus provided 
a more rigorous justification of the Toussaint-Wilczek result for both the t _d//4 density 
decay and amplitude for 2 < d < 4. For d < d c = 2, a full RG calculation becomes 
necessary. Remarkably, comparisons with exact results for the decay exponent in one 
dimension (TBI [TZl EH] show that this qualitative change in the system does not lead to 
any modification in the form of the asymptotic density decay exponent at d = 2 (and so 
very unlike the case of unequal initial denities). However, actually demonstrating this 
using field theory methods has not yet been accomplished, since this would involve a 
very difficult non-perturbative sum over the initial 'surface' terms. 

Lastly, we also mention related work by Sasamoto and coworkers [HE] where the 
mA + nB — > reaction was studied using field-theoretic techniques, by methods similar 
to those of Ref . [21] • These authors also found a t~ d ^ decay rate independent of m and 
n (provided both are nonzero), valid for d < A/(m + n — 1). 

• The segregated A + B — > reaction, reaction zones: 

Two-species reactions can also be studied starting from an initial condition of a 
segregated state, where a (d — l)-dimensional surface separates the two species at time 
t = 0. Later, as the particles have an opportunity to diffuse into the interface, a reaction 
zone forms. Galfi and Racz first studied these reaction zones within the local mean-field 
equations, and were able to extract some rich scaling behavior: the width of the reaction 
region grows asw~ t 1//6 , the width of the depletion region grows, as might be expected, 
as t 1 / 2 , and the particle densities in the reaction zone scale as t" 1 / 3 [H]. 

Redner and Ben-Nairn [EH] proposed a variation of this model where equal and 
opposite currents of A and B particles are directed towards one another and a steady- 
state reaction zone is formed. In this case it is of interest to study how the various 
lengths scale with the particle current J. Within the local mean-field equations they 
found that the width of the reaction region grows as w ~ J -1 / 3 , whereas the particle 
densities in the reaction zone scale as J 2 / 3 . The above initially segregated system may be 
directly related to this steady-state case by observing that, in the former, the depletion 
region is asymptotically much larger than the reaction zone itself. This means there is a 
significant region where the density evolves only by diffusion, and goes from a constant 
to zero over a range L ~ t 1 / 2 . Since J ~ — Va, we find J ~ t _1//2 , which may be used 
to translate results between these two cases. 

Cornell and Droz jHO] extended the analysis of the steady-state problem beyond the 
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mean field equations and, with RG motivated arguments, conjectured a reaction zone 
width w ~ J -1 /^ 1 ) f or th e case d < 2. Lee and Cardy confirmed this result using 
RG methods [01]. The essential physics here is that the only dimensional parameters 
entering the problem are the reaction rate and the current J. However, for d < 2, RG 
methods demonstrate that the asymptotics are independent of the reaction rate. In that 
case, dimensional analysis fixes the above scaling form (with logarithmic corrections 
in d = 2 [21])- Howard and Cardy [02] provided explicit calculations for the scaling 
functions. However, numerical investigations of the exponent of the reaction zone width 
revealed a surprisingly slow convergence to its predicted value w ~ J -1 / 2 ~ t 1 / 4 in 
d = 1 [93J. The resolution of this issue was provided in Ref. [05; where the noise- 
induced wandering of the front was considered (in contrast to the intrinsic front profile 
analyzed previously). There it was shown that this noise-induced wandering dominates 
over the intrinsic front width and generates a multiplicative logarithmic correction to 
the basic w ~ J -1 / 2 ~ t 1 / 4 scaling in d = 1. 

Remarkably, one can also study the reaction zones in the initially mixed system with 
equal initial densities, since it asymptotically segregates for d < 4 and spontaneously 
forms reaction zones. As shown by Lee and Cardy [01], if one assumes that in the 
depletion regions, where only diffusion occurs, the density goes from the bulk value 
t~ d / 4 to zero in a distance of order t 1//2 , the current scales as J ~ £^( d + 2 )/ 4 . From this 
the scaling of the reaction zone width with time immediately follows. As d —>■ 4 from 
below, the reaction zone width approaches t 1//2 , i.e., the reaction zone size becomes 
comparable to the depletion zone, consistent with the breakdown of segregation. This 
analysis also reveals the true critical dimension d c = 2, with logarithmic corrections 
w ~ (tint) 1 / 3 arising from the marginal coupling in d — 2 |21j . 

• Inhomogeneous reactions, shear flow and disorder: 

One important variant of the A+B — > reaction, first analyzed by Howard and Barkema 
[OS], concerns its behavior in the linear shear flow v = t>o?/x, where x is a unit vector 
in the x-direction. Since the shear flow tends to enhance the mixing of the reactants, 
we expect that the reaction kinetics will differ from the homogeneous case. A simple 
generalization of the qualitative arguments of Toussaint and Wilczek shows that this is 
indeed the case. The presence of the shear flow means that in volumes smaller than 
(Dt) d / 2 [l + (wot) 2 /3] 1//2 only the species which was initially in the majority of that region 
will remain. Hence, we immediately identify a crossover time t c ~ Vq . For t <C t c the 
shear flow is unimportant and the usual t~ d l A density decay is preserved. However 
for t t c , we find a t~( d + 2 )/ 4 decay holding in d < 2. Since d = 2 is clearly the 
lowest possible dimension for such a shear flow, we see that the shear has essentially 
eliminated the non-classical kinetics. These arguments can be put on a more concrete 
basis by a field-theoretic RG analysis [OH] , which shows the shear flow adds terms of the 
form avoyd x a and hv^ydj) to the action. The effect of these contributions can then be 
incorporated into modified propagators, after which the analysis proceeds similarly to 
the homogeneous case |2T] . 



Renormalization Group Methods for Reaction- Diffusion Problems 



37 



The related, but somewhat more complex example of A + B — > in a quenched 
random velocity field was considered by Oerding [96 . In this case it was assumed that 
the velocity at every point r = (x, y) of a <i-dimensional system was either parallel or 
antiparallel to the x axis and depended only on the coordinate perpendicular to the flow. 
The velocity field was modeled by quenched Gaussian random variables with zero mean, 
but with correlator (v(y)v(y')) = fo5(y — y'). In this situation, qualitative arguments 
again determine the density decay exponent. Below three dimensions, a random walk 
in this random velocity field shows superdiffusive behavior in the x-direction jHZ|. The 
mean-square displacement in the x direction averaged over configurations of v(y) is 
(x 2 ) ~ £( 5_d )/ 2 for d < 3. Generalizing the Toussaint-Wilczek argument then gives an 
asymptotic density decay of t~^ d+3 ^ 8 . In this case, the system still segregates into A and 
B rich regions, albeit with a modified decay exponent for d < 3. However, to proceed 
beyond this result, Oerding applied RG methods to confirm the decay exponent and 
also to compute the amplitude of the density decay to first order in e = 3 — d (HHl • The 
analysis proceeds along the same lines as the homogeneous case , particularly in the 
derivation of effective 'initial' interaction terms, although care must also be taken to 
incorporate the effects of the random velocity field, which include a renormalization of 
the diffusion constant. Lastly, we mention work by Deem and Park, who analyzed the 
properties of the A + B — > reaction using RG methods in the case of long-ranged 
potential disorder [OH], and in a model of turbulent flow [TTIj . 

• Reversible reactions, approach to equilibrium: 

Rey and Cardy jM] studied the reversible reaction-diffusion systems A + A ^ C 
and A + B # C using RG techniques. Unlike the case of critical dynamics in 
equilibrium systems, the authors found that no new nontrivial exponents were involved. 
By exploiting the existence of conserved quantities in the dynamics, they found that, 
starting from random initial conditions, the approach of the C species to its equilibrium 
density takes the form At~ d l 2 in both cases and in all dimensions. The exponent follows 
directly from the conservation laws and is universal, whereas the amplitude A turns out 
to be model-dependent. Rey and Cardy also considered the cases of correlated initial 
conditions and unequal diffusion constants, which exhibit more complicated behavior, 
including a nonmonotonic approach to equilibrium. 

5.3. Coupled reactions without active phase 

The mixed reaction-diffusion system A + A— > 0, A + B 0, B + B ^ was first studied 
using field-theoretic RG methods by Howard , motivated by the study of persistence 
probabilities (see section [O} . The renormalization of the theory proceeds again similar 
to the case of A + A — > 0: only the reaction rates need to be renormalized, and this 
can be performed to all orders in perturbation theory. For d < d c = 2, perturbative 
calculations for the density decay rates were only possible in the limit where the density 
of one species was very much greater than that of the other. The density decay exponent 
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of the majority species then follows the standard pure annihilation kinetics, whereas 
the minority species decay exponent was computed to 0(e = 2 — d) 57 . This one- 
loop exponent turned out to be a complicated function of the ratio of the A and B 
species diffusion constants. The calculation of this exponent using RG methods has 
been confirmed and also slightly generalized in Ref. jlOOj . The above mixed reaction- 
diffusion system also provides a good testing ground in which to compare RG methods 
with the Smoluchowski approximation, which had earlier been applied to the same multi- 
species reaction-diffusion system . This is a revealing comparison as the value of the 
minority species decay exponent is non-trivial for d < 2, and is no longer fixed purely 
by dimensional analysis (as is the case for the pure annihilation exponent for d < 2). 
This difference follows from the existence of an additional dimensionless parameter in 
the multi-species problem, namely the ratio of diffusion constants. Nevertheless, in this 
case, it turns out that the Smoluchowski approximation decay exponent is identical to 
the RG-improved tree level result, and provides rather a good approximation in d = 1 

IH7] . However, this is not always the case for other similar multi-species reaction- 
diffusion models, where the Smoluchowski approximation can become quite inaccurate 
(see Refs. |5T | 1100] for more details). 

The same system but with equal diffusion constants was also analyzed using RG 
methods in Ref. |102j . as a model for a steric reaction-diffusion system. As pointed 
out in Refs. [HZ1 I103j . this model has the interesting property that at large times for 
d < d c = 2, the densities of both species always decay at the same rate, contrary to the 
predictions of mean-field theory. This result follows from the indistinguishability of the 
two species at large times: below the upper critical dimension, the reaction rates run to 
identical fixed points. Since the diffusion constants are also equal there is then no way 
to asymptotically distinguish between the two species, whose densities must therefore 
decay at the same rate. The same set of reactions, with equal diffusion constants, was 
used to study the application of Bogolyubov's theory of weakly nonideal Bose gases to 
react ion- diffusion systems |104j . 

Related models were studied in the context of the mass distribution of systems of 
aggregating and diffusing particles |105t I56j. In the appropriate limit, the system of 
Ref. jHEj reduced to the reactions A + A — > A and A + B — > 0. Progress could then be 
made in computing to 0{e = 2 — d) the form of the large-time average mass distribution, 
for small masses. Comparisons were also made to Smoluchowski-type approximations, 
which failed to capture an important feature of the distribution, namely its peculiar form 
at small masses, referred to by the authors of Ref. [SB] as the Kang-Redner anomaly. This 
failure could be traced back to an anomalous dimension of the initial mass distribution, 
a feature which, as discussed in section 12 A\ cannot be picked up by Smoluchowski-type 
approximations. Howard and Tauber investigated the mixed annihilation / 'scattering' 
reactions A + A 0, A + A ^ B + B, B + B ^ A + A, and 5 + 5^0 [23 . In 
this case, for d < 2, to all orders in perturbation theory, the system reduces to the 
single-species annihilation case. Physically this is again due to the re-entrance property 
of random walks: as soon as two particles of the same species approach each other, 
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they will rapidly annihilate regardless of the competing 'scattering' processes, which 
only produce particle pairs in close proximity and therefore with a large probability of 
immediate subsequent annihilation. 

Finally, we mention the multi-species pair annihilation reactions A± + Aj — > 
with 1 < i < j < q, first studied by Ben-Avraham and Redner 106, and more 
recently by Deloubriere and coworkers |1U7| 11081 ESI- For unequal initial densities 
or different reaction rates between the species, one generically expects the same scaling 
as for A + B — > asymptotically (when only the two most numerous, or least reactive, 
species remain). An interesting special case therefore emerges when all rates and initial 
densities coincide. For any q > 2 and in dimensions d > 2 it was argued that particle 
species segregation cannot occur, and hence that the asymptotic density decay rate for 
equal initial densities and annihilation rates should be the same as for the single-species 
reaction A + A — > 0. In one dimension, however, particle segregation does take place 
for all q < oo, and leads to a g-dependent power law ~ t~\<i- l )/ 2 <i for the total density 
|107| I108| IllOj . For q = 2, this recovers the two-species decay ~ i~ 1//4 , whereas the 
single-species behavior ~ t -1 / 2 ensues in the limit q — > oo (since the probability that a 
given particle belongs to a given species vanishes in this limit, any species distinction 
indeed becomes meaningless). Other special situations arise when the reaction rates are 
chosen such that certain subsets of the Ai are equivalent under a symmetry operation. 
One may construct scenarios where segregation occurs in dimensions d > 2 despite the 
absence of any microscopic conservation law [109 . 

A variation on this model has a finite number of walkers N of each species Ai, 
initially distributed within a finite range of the origin. Attention is focused on the 
asymptotic decay of the probability that no reactions have occured up to time t. The 
case of Ni — 1 for all i reduces to Fisher's vicious walkers |lllj . and the case N\ = 1 
and N 2 = n reduces to Krapivsky and Redner's lion-lamb model |112j . Applying RG 
methods to the general case, including unequal diffusion constants for the different 
species, Cardy and Katori demonstrated that the probability decays as t~ a " Ni " for 
d < 2, and calculated the exponent to second order in an e = 2 — d expansion |113j . 

5-4- Persistence 

Persistence, in its simplest form, refers to the probability that a particular event has 
never occurred in the entire history of an evolving statistical system |114j . Persistence 
probabilities are often universal and have been found to be nontrivial even in otherwise 
well- understood systems. An intensively studied example concerns the zero-temperature 
relaxational dynamics of the Ising model, where one is interested in the persistence 
probability that, starting from random initial conditions, a given site has never been 
visited by a domain wall. In one dimension, the motion and annihilation of Ising domain 
walls at zero temperature is equivalent to an A + A — > reaction-diffusion system, where 
the domain walls in the Ising system correspond to the reacting particles. An exact 
solution exists for the persistence probability in this case |115j . but, as usual, the solution 



Renormalization Group Methods for Reaction- Diffusion Problems 



40 



casts little light on the question of universality. A different approach was proposed by 
Cardy who studied, in the framework of the reaction-diffusion model, the proportion of 
sites never visited by any particle In d = 1 (though not in higher dimensions) this 

is the same quantity as the original persistence probability. Furthermore, since Cardy 
was able to employ the kind of field-theoretic RG methods discussed in this review, the 
issue of universality could be addressed as well. 

Cardy demonstrated that the probability of never finding a particle at the origin 
could be calculated within the field-theoretic formalism through the inclusion of an 
operator product fit <5 a t ao Q . The subscript denotes that the dld operators are associated 
with the origin, and the operator-valued Kronecker 5-function ensures that zero weight 
is assigned to any histories with a particle at the origin. This operator has the net effect 
of adding a term —h Jq 0(0, t')(j)(0, t')dt' to the action, and the persistence probability 
then corresponds to the expectation value (exp(— h Jq 0(0, t')dt')), averaged with respect 
to the modified action. Power counting reveals that [h] ~ K, 2 ~ d , so this coupling is 
relevant for d < 2. Cardy showed that renormalization of this interaction required both 
a renormalized coupling and a multiplicative renormalization of the field 0(0, t). 
This results in a controlled e = 2 — d expansion for the universal persistence exponent 
6 = 1/2 + 0(e) |llb'j . This compares to the exact result in one dimension by Derrida et 
al., namely 8 = 3/8 |115j . An alternative approach to this problem was given by Howard 
(HZ] in the mixed two-species reaction A + A ^ 0, A + B 0, with immobile B particles 
(see also section UTS)) . In this case the persistence probability corresponds to the density 
decay of immobile B particles in d = 1, in the limit where their density is much smaller 
than those of the A particles. Howard's expansion confirmed the results of Cardy and 
also extended the computation of the persistence probability to 0(e = 2 — d). The case 
of persistence in a system of random walkers which either coagulate, with probability 
(q — 2)/(q — 1), or annihilate, with probability l/(q — 1), when they meet was also 
investigated using RG methods by Krishnamurthy et al. |117j . In one dimension, this 
system models the zero-temperature Glauber dynamics of domain walls in the g-state 
Potts model. Krishnamurthy et al. were able to compute the probability that a given 
particle has never encountered another up to order e = 2 — d. 

A further application of field-theoretic methods to persistence probabilities was 
introduced by Howard and Godreche |f 18j in their treatment of persistence in the voter 
model. The dynamics of the voter model consist of choosing a site at random between 
t and t + dt; the 'voter' on that site, which can have any of q possible 'opinions', then 
takes the opinion of one its 2d neighbours, also chosen at random. This model in 
d = 1 is identical to the Glauber-Potts model at zero temperature, but can also, in 
all dimensions, be analyzed using a system of coalescing random walkers. This again 
opens up the possiblity for field-theoretic RG calculations, as performed in Ref. |118j . 
The persistence probability that a given 'voter' has never changed its opinion up to 
time t was computed for all d > 2, yielding an unusual exp[— /(g)(lnt) 2 ] decay in two 
dimensions. This result confirmed earlier numerical work by Ben-Nairn et al. |119j . 
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6. Active to Absorbing State Transitions 

In the previous sections, we have focused on the non-trivial algebraic decay towards the 
absorbing state in diffusion-limited reactions of the type kA — > iA (with k > 2 and 
£ < k), and some variants thereof. Universal behavior naturally emerges also near a 
continuous nonequilibrium phase transition that separates an active state, with non- 
vanishing particle density as t — > oo, from an inactive, absorbing state. We shall see 
that generically, such phase transitions are governed by the power laws of the directed 
percolation (DP) universality class [1271111211 l£l 112] . 



6.1. The directed percolation (DP) universality class 

A phase transition separating active from inactive states is readily found when 
spontaneous particle decay (A — > 0, with rate ji) competes with the production 
process (A — > A + A, branching rate cr). In this linear reaction system, a(t) = 
a(0) exp[— (/i — a)t] — > exponentially if a < fi. In order to render the particle density 
a finite in the active state, i.e., for a > fi, we need to either restrict the particle number 
per lattice site (say, to or 1), or add a binary reaction A + A — > (0, A), with rates 
A(A'). The corresponding mean- field rate equation reads 

d t a(t) = (a - fi)a(t) - (2A + A') a(t) 2 , (68) 
which for a > fi implies that asymptotically 

a{t) ^a^= 2 ^ ^, , (69) 

which is approached exponentially \a(t) — a^l ~ exp[— (cr — n)t] as t — > oo. Precisely at 
the transition a = fi, Eq. ()68|) yields the binary annihilation / coagulation mean-field 
power law decay a(t) ~ t -1 . Generalizing Eq. (jfi8j) to a local particle density and taking 
into account diffusive propagation, we obtain with r = (yU — o~)/D: 

d t a(x, t) = -D(r - V 2 )a(x, t) - (2A + A') a(x, tf , (70) 

wherefrom we infer the characteristic length and diffusive time scales £ ~ |r| -1 / 2 and 
t c ~ C, 2 /D ~ which both diverge upon approaching the critical point at r = 0. 
Upon defining the critical exponents 

(<0~(-rf (r<0), (a(t))~r a (r = 0) , 

e~|r|^ (r^0), t c ~C/D~\r\- z » (r ^ 0) , (71) 

we identify the mean-field values (3 — 1, a — 1, u — 1/2, and z = 2. 

In order to properly account for fluctuations near the transition, we apply the field 
theory mapping explained in sectional The ensuing coherent-state path integral action 
then reads 

S[4>, 0] = J d d x{-<p{t f ) + J* f dt [4> (d t - DV 2 ) <f>-n{l- 4>)<j) + a(l - 0)# 

-A (l - 2 ) 2 - A' (l - 0) 00 2 ] - n o 0(O)|, (72) 
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which constitutes a microscopic representation of the stochastic processes in question. 
Equivalently, we may consider the shifted action (with = 1 + 0) 

5[0, 0] = J d d x J dt {0 [d t + D(r - V 2 )] - (J0 2 + (2A + A')00 2 + (A + A')0 2 2 } . (73) 

Since the ongoing particle production and decay processes should quickly obliterate any 
remnants from the initial state, we have dropped the term no0(O), and extended the 
temporal integral from — oo to oo. The classical field equations 5S/5<ft = (always 
solved by = 0) and 5S/5(j) = yield the mean- field equation of motion (fTUj) . 

Our goal is to construct an appropriate mesoscopic field theory that captures the 
universal properties at the phase transition. Recall that the continuum limit is not 
unique: We are at liberty to choose the scaling dimensions of the fluctuating fields 
0(x, t) and 0(x, t), provided we maintain that their product scales as a density, i.e., 
[00] = K d with arbitrary momentum scale k. In RG terms, there exists a redundant 
parameter |122j that needs to be eliminated through suitable rescaling. To this end, we 
note that the scaling properties are encoded in the propagator G(x, t) = (0(x, t)0(O, 0)). 
The lowest-order fluctuation correction to the tree-level expression 

G (p,u>)= . * - 2 , (74) 
— iuj + D[r + p z ) 

is given by the Feynman graph depicted in figure Efb, top), which involves the product 
~ — cr(2A + A') of the two three-point vertices in (|7Hjl . Similarly, the one-loop correction 
to either of these vertices comes with the very same factor. It is thus convenient 
to choose the scaling dimensions of the fields in such a manner that the three-point 
vertices attain identical scaling dimensions. This is achieved via introducing new fields 



s(x,t) = 0(x,tW(2A + X')/cr and s(x,i) = 0(x, t)^a/(2X + A'), whence 

S[s, s}= J d d x J dt [s [d t + D(r - V 2 )] s - u(s - s)ss + (A + A')s 2 s 2 } . (75) 

Here, u = \Jo~{2\ + A') is the new effective coupling. Since [a] = k 2 and [A] = n 2 ~ d = 
[A'], its scaling dimension is [u] = n 2 ~ d / 2 , and we therefore expect d c — 4 to be the upper 
critical dimension. Moreover, [(A + X')/u] = K~ d l 2 scales to zero under subsequent RG 
transformations: compared to u, both couplings A and A' alone constitute irrelevant 
parameters which will not affect the leading universal scaling properties. 

Upon omitting these irrelevant terms, we finally arrive at the desired effective field 
theory action 

S cS [s, s} = J d d x J dt [s [d t + D(r - V 2 )] s - u(s - s)ss} . (76) 

It displays duality invariance with respect to time (rapidity) inversion, s(x, t) <-> 
— s(x, —t). Remarkably, the action ()76|) was first encountered and analyzed in particle 
physics under the guise of Reggeon field theory \l'23[ 1124] . It was subsequently noticed 
that it actually represents a stochastic ('Gribov') process |125tll26j . and its equivalence 
to the geometric problem of directed percolation was established [12*3 I128[ I120j . In 
directed bond percolation, randomly placed bonds connecting regular lattice sites can 
only be traversed in a given preferred special direction, which is to be identified with t 
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in the dynamical problem. Particle decay, coagulation, and production respectively 
correspond to dead ends, merging links, or branching of the ensuing percolating 
structures. Near the percolation threshold, the scaling properties of the critical 
percolation cluster are characterized by the exponents governing the divergences of 
the transverse correlation length u± — v and of the longitudinal (in the t direction) 
correlation length u» = zv. (For more details, we refer the reader to Refs. |Hl I12j.) 

From our derivation of the effective action (|75jl above, it is already apparent that 
either pair annihilation or coagulation lead to identical critical properties. Instead of 
these binary reactions, we could also have employed site oocupation number restrictions 
to render the particle density finite in the active phase. Van Wijland has recently 
shown how such local constraints limiting to values of 0, 1 only can be implemented 
into the second-quantized bosonic formalism |129j . thus avoiding a more cumbersome 
representation in terms of spin operators. The resulting action acquires exponential 
terms for each field (p. For the competing first-order processes A — > (0, 2 A) one 
eventually obtains 

5 rest [0,0] = J d d x J rft[-Ml-0)0e-^ + a(l-0)#e- 2 ^] , (77) 

where we have merely written down the bulk reaction part of the action, and v is a 
parameter of scaling dimension [v] = n~ d which originates from taking the continuum 
limit. Since therefore v will scale to zero under RG transformations, we may expand the 
exponentials, whereupon the leading terms in the corresponding shifted action assume 
the form (jTKjl. with 2A + A' = {2a — fi)v ~ ov and A + A' = 4crt>. Thus we are again led 
to the effective DP field theory action ()76|) (despite the formally negative value for A). 

Following the procedure outlined in section 13*31 I], we find that the field theory 
action (j7fij) is equivalent to the stochastic differential equation 

d t s = D (v 2 -r) s -us 2 + V2u~sr] , (78) 

with (77) = 0, (77 (x, t) 77 (x', t')) = 5(x — x') 6(t — t'), or, upon setting ( = \/2usi] 
in order to eliminate the square-root multiplicative noise, (() = 0, (£(x, i) £(x', t')) = 
2ws(x, £)5(x— x') 5(t—t f ). We may view these resulting terms as representing the leading- 
order contributions in a power-law expansion of the reaction and noise correlation 
functionals R[s] = r + us + . . . and N[s] = u + . . . with respect to the density s of 
activity in 

d t s = D (V 2 — R[s}) s + C, (C(x, t) C(x , 0) = 2sN[s}5(x - x') 5(t - t') , (79) 

which represents the general Langevin description of systems displaying active and 
absorbing states ^2] • A factor s has been factored out of both R and N here, since the 
stochastic processes must all cease in the inactive, absorbing phase. These considerations 
establish the DP hypothesis: The critical properties near an active to absorbing state 
phase transition should generically be governed by the directed percolation scaling 
exponents, provided the stochastic process is Markovian, the order parameter decoupled 
from any other slow variable, there is no quenched disorder in the rates, and no special 
symmetries require that any of the lowest-order expansion coefficients r or u vanish 
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|12U| I121| . There is even a suggestion that the glass transition in supercooled liquids 
might be governed by a zero-temperature fixed point, with critical exponents in the DP 
universality class 




Figure 7. DP field theory: (a) vertices, and (b) one-loop Feynman graphs for the 
two- and three-point vertex functions. 



6.2. Renormalization and DP critical exponents 

The asymptotic scaling behavior of DP can be inferred from the renormalized propagator 
G(p, (J) = r^^p, — uj)~ 1 . The tree contribution is given by ([74ft : by combining the two 
three-point vertices in figure Efa) one arrives at the one-loop Feynman graph depicted 
in figure [T^b, top), whose corresponding analytic expression reads in Fourier space 

„ 2 f d d p' rduj' 1 
2u ' 1 



{2Ti) d J 2vr -i(u' + u/2) + D[r + (p' + p/2) 2 ] 

1 

_,(_^ + u / 2 ) + D[r + (-p' + p/2) 2 ] ' (80) 

if we split the external momentum and frequency symmetrically inside the loop. The 
integration over the internal frequency u' is now readily performed by means of the 
residue theorem, whereupon we obtain 

.2 , d d p , 1 

(2vr) d iu/2D + r + p 2 /A + p' 2 
The loop contribution displays IR singularities as r — > 0, u — > 0, and p — ► 0. In the 
ultraviolet, it diverges in dimensions d > 2. The leading divergence, however, can be 
absorbed into a fluctuation-induced shift of the critical point away from the mean-field 
r = 0. On physical grounds one must demand G(p = 0,u = 0) _1 = at criticality. 
Consequently, the new critical point is given self-consistently by 

u 2 r d d p' 1 



2 

r( 1 ' 1 )(p,^) = ^ + J D(r + J 9 2 ) + ^ J 



+ 0(u 4 ) . (82) 



D 2 J (2vr) d r c + p' 2 

Fluctuations tend to increase the likelihood of extinction (if the density is already 
low, a chance fluctuation may drive the system into the absorbing state), and thus 
reduce the parameter regime of the active phase as compared with mean-field theory. 
In dimensional regularization, one assigns the value 

I(r\=I**- 1 = Vi ~ S ~ d/2) r-^l 2 (83) 
s[ ' J (2n) d (r+p 2 Y 2 d Ti d / 2 T(s) ' 1 ' 
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also to those momentum integrals that are UV- divergent. The solution to (J82|) then 
reads explicitly \r c \ = [2A d u 2 /(d - 2) (4 - d)D 2 ] 2 /^ with A d = T(3 - d/2)/2 d - 1 n d / 2 . 
The shift of the transition point thus depends nonanalytically on e = 4 — d. 

Let us introduce the true distance from the critical point r = r — r c . Upon inserting 
JB2} into (JHIJ), we find 

. 9x w 2 /• dV iu/2D + r + p 2 /4 ^. 4x . 

DJ (2ix) d p> 2 Uu/2D + t + p 2 /4 + p l2 \ 

The integral here is UV-divergent in dimensions d > 4. There are three such singular 
terms, proportional to iu, Dt, and -Dp 2 , respectively. Consequently, we require three 
independent multiplicative renormalization factors to render the two-point function or 
propagator finite. In addition, the three-point vertex functions I^ 1,2 ) and r( 2,1 ) carry 
(identical) UV-singularities for d > 4. We thus define renormalized parameters D R , t r , 
and u R , as well as renormalized fields s R according to 

s R = Z]/ 2 s , D R = Z D D , t r = Z T r k- 2 , u R = Z u uA 1 ' 2 k^' 2 . (85) 

As a consequence of rapidity inversion invariance, s R = Z x J 2 s as well, whence T R ^ = 
Z~ 1 T( 1 ' 1 \ In the minimal subtraction prescription, the Z factors contain merely the 1/e 
poles with their residues. Choosing the normalization point r R — 1, u — 0, p = 0, we 
may read off Z s and the products Z s ZdZ t , Z s Zd from the three terms on the right-hand 
side of (jEU), and therefrom to one- loop order 

u 2 A d K~ € u 2 A d K~ € 3u 2 A d K~ e , , 

This leaves just Z u to be determined. It is readily computed from the three-point 
function r'- 1 ' 2 -', whose one-loop graph is depicted in figure Efb, bottom), or from I^ 2 ' 1 ). 
At the normalization point (NP), 

r <M,| NP = _ r <M)| NP = _ 2u U - % I** —Lj 

\ D 2 J (2vr) d (r + p 2 ) , 



. 2u(l _^^L], (87) 



which directly yields the product Z^ 2 Z U , and with 

5u 2 A d nT e 



(88) 



AD 2 e 

Since all higher vertex functions are UV-finite, this completes the renormalization 
procedure for the DP field theory (J76j) . We identify the effective coupling constant 
as v = u 2 /D 2 , with renormalized counterpart 

v R = Z v vA d K d ~ A , Z v = Z 2 JZ 2 D . (89) 

We may now write down the DP generalization of the Callan-Symanzik equation 
(foTIJ) for the propagator, recalling that G R = Z S G: 



d Odd 
Cs + CdD r — — + CtTr -= \- Pv (vr) 



ok dD R dr R dv R 



Gr(P, U, D r , Tr, K, Vr) = , (90) 
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(91) 
(92) 
(93) 
(94) 



In dimensions d < d c — 4 (e > 0), the (3 function yields a nontrivial stable fixed 
point 



^=3+0(6 2 ) 



(95) 



Solving ()90j) with the method of characteristics, /t — > and using the form (|74jl . we 
find in its vicinity the scaling law 



G R (p, u, D R , tr, k, vr)' 1 ~ p 2 £) fl ^.(«y+fc(*i) f 



to 



D r £^ v *r\k£) 2 



,tr£^\v r ,(96) 



with T representing a dimensionless scaling function. Upon employing the matching 
condition £ 2 = p 2 /k 2 , this yields 



Gr(p,u,D r ,t r ,k,v* r 



D^ 2 |p| 2 -"f 1, 



^IpI 



:,Tr P 



with the t/iree independent scaling exponents 

V = -Cs(v R )-CD(v R ) = ~^ + 0(e 2 ) 

z = 2 + ( D (v R ) = 2-^ + 0(e 2 ) , 

= -Uv* R ) = 2-- + 0(e 2 ) . 



7/ 



4 



Alternatively, with £ = \r R \ u we arrive at G# ~ |r| 7 , where 

7 = i/(2 - ,7) = 1 + | + Ofe 3 ) . 

o 



v*r ) , (97) 

(98) 
(99) 
(100) 

(101) 



In a similar manner, we obtain in the active phase for the average {s R ): 

(s R {t,D R ,r R ,K,v R )) = , (102) 



d C s d d d 

K 7T + (dD r — — + ( t t r — + (3 v {v R ) — 

OK 2 ODr OTr OVr 



whose solution near the stable fixed point v* R reads 

(s R (t,D R ,T R ,K 7 v R )) ~ K^l^W ^Dntfe^^TBpW^ 

where (s R ) = K d l 2 s. Consequently, matching i 
respectively, gives (s R ) ~ |trP and (s R ) ~ t~ a , with 



\t r \ u and 



P 



v[d-(.(v* R )] v(d + ri + z-2) 



P 



l-^ + 0(e 2 ) 



a = — = 1 - - + 0(e 2 ) 
zv A y ' 



(103) 

(104) 
(105) 
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Field-theoretic tools in conjunction with the RG therefore allow us to define the 
generic universality class for active to absorbing state phase transitions, derive the 
asymptotic scaling laws in the vicinity of the critical point, and compute the critical 
exponents perturbationally by means of a systematic expansion about the upper critical 
dimension d c = 4. In higher dimensions d > 4, the only fixed point is Vr = 0, and 
we recover the mean-field scaling behavior. Precisely at the upper critical dimension, 
there appear logarithmic corrections. These, as well as the two-loop results for the 
critical exponents and the scaling behavior of various other observables, are presented 
in Ref. Monte Carlo simulations have determined the numerical values for the 

DP critical exponents in dimensions d < 4 to high precision [HIE], and confirmed the 
logarithmic corrections predicted by the RG jl31|. I132j . 

6. 3. Variants of directed percolation processes 
• Multi-species DP processes: 

We argued in section 16.11 that absorbing to active phase transitions should generically 
be described by the critical exponents of DP. This far-reaching assertion is based on 
the structure of Eq. (f7T?|) . identifying the field s as some coarse-grained 'activity' density 
|120|I121] . It is indeed very remarkable that the DP universality class extends to multiple 
species of reacting agents. Consider, for example, the reactions A # A + A, A — > 0, 
coupled to a similar system B ^ B + B, B — > via the processes A — > B + B, 
A + A —>■ B, and its obvious extension to additional reactants. Inclusion of higher- 
order reactions turns out not to change the critical properties, since the corresponding 
couplings are all irrelevant under the RG. One then arrives at the following effective 
Langevin description for coupled coarse-grained density fields [1331 EH] : 

dtst = A (V 2 - Ri [ Si ]) Si + Q, Ri [ Si ] = n + £ g i3 s 3 + ... (106) 

j 

(Ci(x,t) 0(x',t')> = 2s i N i [s i )8 ij 6{x - x')5(t - t') , Ni[ Si ] =m + ..., (107) 

generalizing Eq. (|79|). As Janssen has demonstrated, the ensuing renormalization factors 
are all given precisely by those of the single-species process, whence the critical point is 
generically described by the ordinary DP scaling exponents. 

Yet if first-order particle transmutations A — > B, etc. are added (notice these 
are also effectively generated by the above reactions), leading to additional terms 
~ J2j^i9jSj in Eq. ()106j) . one finds that the ensuing RG flow typically produces 
asymptotically unidirectional processes. One then encounters multicritical behavior, 
if several control parameters vanish simultaneously |134[ 1135] 11361 ITTTj . While the 
critical exponents rj, v, z, and 7 remain unchanged, there emerges in this situation 
a hierarchy of order parameter exponents (3k — l/2 fc — 0(e) on the kth level, e.g., 
p 1 = (3 = 1 - e /6 + 0(e 2 ), fa = 1/2 - 13e/96 + 0(e 2 ), and similarly for the decay 
exponents = ftk/zv. The crossover exponent associated with the multicritical point 
can be shown to be = 1 to all orders in the perturbation expansion [HH] • 
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• Dynamic isotropic percolation (dIP): 

An alternative mechanism to generate novel critical behavior in a two-species system 
operates via a passive, spatially fixed and initially homogeneously distributed species 
X that couples to the diffusing and reproducing agents A — > A + A through the decay 
processes A — > X and X + A — > X. Upon integrating out the fluctuations of the inert 
species X, and expanding about the mean-field solution, the resulting effective action 
eventually becomes 

S e «[s, s} = J d d x J dt Ss [d t + D(r - V 2 )] s - us 2 s + Duss J* s(t')dt'\ , (108) 
which corresponds to a stochastic differential equation 

d t s = D(V 2 -r)s- Dus [ s(t')dt' + V2u~sr] , (109) 

^ ' J —oo 

with (77) = 0, (r){x,t)T]('x!,t')) = <5(x — x') 5(t — t') as in DP. Thus, the induced decay 
rate is proportional to the product of the densities of active agents s and the 'debris' 
D J^oo s(t')dt' produced by decayed agents A. More generally, the action ()108|) describes 
the general epidemic process \1'SS\ 1140} IT4*T] . in contrast with the simple epidemic 
process represented by DP. 

We may now consider the quasistatic limit of the field theory ()108J) by introducing 
the fields (p — s(t — > 00), and = D Jf^ s(t')dt', whereupon we arrive at the action 

Sqst[<p, <p] = jd d x<p r - V 2 - u(<p - if) ^. (110) 

Note, however, that as a manifestation of its dynamic origin, this quasi-static field 
theory must be supplemented by causality rules. The action (jllOj) then describes the 
scaling properties of critical isotropic percolation clusters |142j . Thus, as first remarked 
by Grassberger |lH8j . the general epidemic process is governed by the static critical 
exponents of isotropic percolation. These are readily obtained by means of a RG analysis 
in an e expansion about the upper critical dimension d c = 6. The one- loop diagrams are 
precisely those of figure [3 with the static propagator Go(p) = ( r +p 2 ) -1 . The explicit 
computation proceeds as outlined in section lo^ (for more details, see, e.g., Ref. [T2]). 
and yields 77 = -e/21 + 0(e 2 ), v' 1 = 2 - 5e/21 + 0(e 2 ), and /3 = 1 - e/7 + 0(e 2 ), with 
e = 6 — d. 

In order to characterize the dynamic critical properties, however, we must return to 
the full action ()108j) . Yet its structure once again leads to the Feynman graphs depicted 
in figure [7[ but with the second vertex in (a) carrying a temporal integration. One then 
finds z = 2 — e/6 + 0(e 2 ), which completes the characterization of this dynamic isotropic 
percolation (dIP) universality class jl39|. I140|. I141j . Precisely as for DP processes, multi- 
species generalizations generically yield the same critical behavior, except at special 
multicritical points, characterized again by a crossover exponent 0=1 



• Levy flight DP: 

Long-range interactions, as can be modeled by Levy flight contributions D^p u to the 
propagators, may modify the critical behavior of both DP and dIP |%2"| Two 
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situations must be distinguished |143j : for 2 — a = 0(e), a double expansion with 

respect to both e and 2 — a is required; on the other hand, if 2 — a = 0(1), the 

ordinary diffusive contribution Dp 2 to the propagator becomes irrelevant, whereas the 

non-analytic term D^p 17 acquires no fluctuation corrections, whence Z s Zd l = 1 to all 

orders in perturbation theory. Subsequent scaling analysis yields the critical dimensions 

d c = 2a for DP and d c = 3a for dIP, respectively. To one-loop order, one then finds at the 

new long-range fixed points, with e = d — d c : rj = 2 — a, and for DP: z = a — e/7 + 0(e 2 ), 
jy-i = (j 0^/7 i r\( A\ a — i _ 0^/7^- i m A\- f™. ^itd. „ — ^ _ q^/i« i /0/V2 



7^ 



-1 



2e/7 + 0(e 2 ), /3 = 1 - 2e/7a + 0(e 2 ); for dIP: z = a - 3e/16 + 0(e 2 
a - e/4 + 0(e 2 ), /3 = 1 - e/4a + 0(e 2 



• DP coupled to a non-critical conserved density (DP-C): 

A variant on the DP reaction scheme A ^ A + A with decay A — > is to require the 
A particle decay to be catalyzed by an additional species C, via A + C — > C. The C 
particles move diffusively and are conserved by the reaction. In the population dynamics 
language, the C particles can be said to poison the A population [144 . Related is a 
model of infection dynamics, A+B — > 2B, B — > A, where A and B respectively represent 
healthy and sick individuals |145[ I14f>j . The latter model reduces to the former in the 
case of equal diffusion constants Da = D b (see Ref. [12 for details). 

These systems exhibit, like DP, an upper critical dimension d c = 4. Below this 
dimension there exist three different regimes, depending on whether the ratio of diffusion 
constants greater than, equal to, or less than unity. For the case Da > D B , the resulting 
RG flows run away, indicating a fluctuation-induced first-order transition |146j . For 
the case Da = D b the critical exponents are given by z = 2, u = 2/d (exact), and 
P = 1 — e/32 + 0(e 2 ); while for Da < D B a distinct fixed point is obtained |145j with 
exact values z — 2, v — 2/d, and [3=1. 

6.4- Branching and annihilating random walks (BARW) 

Branching and annihilating random walks are defined through diffusing particles A 
subject to the competing branching reactions A — > (m + 1)A (with rate a) and 
annihilation processes kA — ► (rate A) |147j . The corresponding mean-field rate 
equation for the particle density reads 

d t a(t) = aa(t) - k\a(t) k , (111) 

with the solution 

«W = — 7777— TT , (112) 

which for t ^> I /a approaches the finite density = (a/kX) 1 /^ 1 ^, independent of the 
initial density ao- Mean-field theory therefore predicts only an active phase, provided 
a > 0. At a c = 0, there exists a 'degenerate' critical point, whose critical exponents are 
given by the pure diffusion-limited annihilation model, r] = 0, v = 1/2, z = 2, 7 = 1, 
a = (3= l/(ife- 1). 
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This mean-field scenario should hold in dimensions d > d c = 2/(k — 1), as 
determined from the scaling dimension of the annihilation process. The branching rate 
with scaling dimension [a] = k 2 constitutes a relevant parameter. However, as initially 
established in Monte Carlo simulations j!47| El Ej, for k — 2 fluctuations invalidate 
this simple picture, rendering BARW considerably more interesting |148[ 1149] . Indeed, 
one has to distinguish the cases of odd and even offspring numbers m: Whereas the 
active to absorbing transition in BARW with odd m is described by the DP critical 
exponents, provided a c > 0, BARW with even m define a genuinely different parity- 
conserving (PC) universality class, named after its unique feature of conserving the 
particle number parity under the involved reactions. 

The field-theoretic representation for BARW with k = 2 and m offspring particles, 
omitting the temporal boundary terms, reads 

S[$, 0] = J d d x J dt [j> (d t - DV 2 ) <p + a(l- 4> m )(f)(j) - A (l - 2 ) (j) 2 ] , (113) 

and the corresponding vertices are shown in figure Efa). Upon combining the branching 
with the pair annihilation processes as in the top one-loop diagram in figure Efb), we 
see that in addition to the original A — > (m + 1) A reaction all lower branching processes 
A — > [m — 1) A, (m — 3) A, . . . become generated. In a first coarse-graining step, all these 
reactions then must be added to the 'microscopic' field theory (jllHj) . Furthermore, upon 
inspecting the renormalization of the branching and annihilation rates by the one-loop 
Feynman graphs depicted in figure Efb), we see that identical loop integrals govern 
the corresponding UV singularities, but the renormalization of the branching process 
with m offspring carries a relative combinatorial factor of m(m + l)/2. Since the loop 
contributions carry a negative sign, the resulting downward shift of the branching rate's 
scaling dimension is lowest for m = 1 and m = 2, respectively. For odd offspring number 
m, the most relevant emerging branching process thus is A — > A + A, but in addition 
the spontaneous decay A — > is generated |148| I149j . Consequently, the effective field 
theory that should describe BARW with odd m becomes dZHJ), and the phase transition is 
predicted to be in the DP universality class. This is true provided the induced particle 
decay rate may overcome the generated or renormalized branching rate with single 
offspring, and thus shift the critical point to o c > 0. Within a perturbational analysis 
with respect to the annihilation rate A, this happens only in low dimensions d < 2 
|148| I149j . In a recent nonperturbative numerical RG study, however, the emergence 
of an inactive phase and DP critical behavior was found in higher dimensions as well 

[na [[Hi- 
lt now becomes apparent why BARW with even offspring number should behave 
qualitatively differently: in this case, spontaneous particle decay processes A — > 
cannot be generated, even on a coarse-grained level. This is related to the fact that 
in the branching processes A — > (m + I) A as well as the pair annihilation A + A —>■ 
the particle number parity remains conserved; correspondingly, there are two distinct 
absorbing states for even-offspring BARW, namely the strictly empty lattice if the initial 
particle number N is even, and a single remaining particle, if N is odd. In the field 
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Figure 8. BARW field theory: (a) branching (top) and annihilation (bottom) 
vertices; (b) one-loop Feynman diagrams generating the A — > (m — 1)A process, and 
renormalizing the branching and annihilation rates, respectively. 



theory ()113j) . this conservation law and symmetry are reflected in the invariance with 
respect to simultaneously taking — > — and — > — 0. This invariance must be 
carefully preserved in any subsequent analysis. Performing the field shift 0=1 + 
masks the discrete inversion symmetry; worse, it becomes lost entirely if afterwards, 
based on mere power counting arguments, only the leading powers in are retained, 
whence one would be erroneously led to the DP effective action. It is therefore safest to 
work with the unshifted action (jllHj) . but adding the generated branching processes 
with m — 2,m — 4, ... offspring particles. As explained before, the most relevant 
branching reaction will be the one with two offspring. Setting m = 2 in the action 
f)l 13j) indeed yields a renormalizable theory, namely the effective action for the PC 
universality class, with the particle production processes with higher offspring numbers 
constituting irrelevant perturbations. 

The bare propagator of this theory is similar to Eq. (J73J), 

G {p,u)= . * n 2 (114) 
— iuj + cr + Dp A 

but contains the branching rate a as a mass term. The branching rate also appears 
in the three-point vertex, figure |Hta,top). Since we need to follow the RG flow of the 
renormalized reaction rates 

o R = Z a( r/D k\ X r = Z x \C d /DK 2 ~ d , (115) 

with Cd = T(2 — d/2)/2 d ~ 1 -K d / 2 , we must set the normalization point either at finite 
external momentum p = 2k, or frequency / Laplace transform variable iuo = s = 2Dk 2 . 
From the one-loop Feynman graphs in figure |Hfb) that respectively describe the 
propagator, branching vertex, and annihilation vertex renormalizations, one then finds 
that the UV singularities, for any value of a, can be absorbed into the Z factors 
^ 3C d X/D C d X/D 

2 - d (k 2 + a/Dy- d / 2 ' A 2 - d (k 2 + a/Dy~ d / 2 ' K ' 

which are functions of both a / D and X/D, as in other crossover theories with relevant 
parameters. With ( a = nd K \n.(o- R / a) and (\ = nd K ln(A^/A), we are thus led to the 
coupled RG flow equations 



Renormalization Group Methods for Reaction- Diffusion Problems 



52 



\ R (£)( x (£)=\ R (£)[d-2 + 



l+O-R 



2-d/2 



'1181 



The effective coupling controlling the RG flows, to one-loop order at least, is 
9r — A/j/ (1 + a R ) 2 ~ d ^ 2 . For a R = 0, that is for the pure pair annihilation model, 
according to Eq. ()118j) g R — > g R = 2 — d, which after trivial rescaling corresponds to the 
annihilation fixed point (JoTlj) . For a > 0, however, we expect cr R (£) — > oo, whereupon 
the RG P function for the coupling g R becomes 



P 9 (g 



R) 



9r 



Ca- 2 



9R 



2 - 



10 - 3d 



-9r 



(119) 



which yields the Gaussian fixed point at g R — and a critical fixed point g* = 
4/(10 — 3c?). Yet since the bare reaction rate corresponding to the pure annihilation 
fixed point is already infinite, see section l4~2l following Eq. (joTjl . we must demand on 
physical grounds that g* < 2 — d, whence we infer that the critical fixed point comes 
into existence only for d < d' c = 4/3. If initially g R < g*, g R (£) — > 0, consistent with 
&r(£) ~ t~ 2 — ► °o as t — > 0. This Gaussian fixed point, characterized by naive scaling 
dimensions, describes the active phase with exponential correlations. On the other hand, 
for g R > g*, and provided that d < d' c , o~ R (^) — > and g R (£) — > 2 — d, which describes an 
inactive phase that in its entirety is governed by the pure annihilation model power laws. 
The phase transition in the PC universality class, which apparently has no counterpart 
in mean-field theory, is thus triggered through fluctuations that drive the branching 
rate irrelevant. In contrast to equilibrium systems, fluctuations here open up a novel 
phase rather than destroying it, and we may view the new borderline dimension d' c as an 
'inverted lower' critical dimension, since the phase transition only exists for d < d' c . The 
phase diagram as function of spatial dimension, within the one-loop approximation, is 
summarized in figure 



i/ g 

3 



ACTIVE 
PHASE 



INACTIVE 
PHASE 



MEAN 
FIELD 



Figure 9. Phase diagram and unstable RG fixed point 1/ g* R for even-offspring BARW 
(PC universality class) as function of dimension d (from Ref. 149 ). 



In order to obtain the asymptotic scaling behavior for the particle density, we write 
down the solution of its RG equation in the vicinity of an RG fixed point g* Rl which 
reads 

a(t, D R , a R , X R , k) = K d £ d a[D R K 2 t 2 t, a R £°^\ X R £^) , (120) 
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since there is no renormalization of the fields or diffusion constant to one-loop order, 
which immediately implies 77 = and z — 2. For d' c < d < 2 the branching rate a R 
plays the role of the critical control parameter tr in DP, and g* R = 2 — d = e is the 
annihilation fixed point. In an e expansion about the upper critical dimension d c = 2, 
we thus obtain the scaling exponents 

v~ x = -Ca{g* R ) = 2 - 3e , z = 2, a = d/2, (3 = dv = zua , (121) 

the latter via matching a R i.^^*^ = 1. Notice that v diverges as e — ► 2/3 or d — > rf' c . 

Yet the PC phase transition at a c > can obviously not be captured by such an e 
expansion. One is instead forced to perform the analysis at fixed dimension, without the 
benefit of a small expansion parameter. Exploiting the mean-field result for the density 
in the active phase, for d < d' c we may write in the vicinity of g*: 

a(t, D R , a R , X R , k) = n d ^ ~ a (a R tn 2 i 2+ ^\ e R ^A , (122) 

where e oc g* — g constitutes the control parameter for the transition, and ( £ = d(3 g /dgR. 
Now setting e R f^ = 1, we obtain with Q a {g* c ) = -2(4 - 3d)/(10 - 3d), ( x (g*) = 
-(4 - d)(4 - 3d)/(10 - 3d), and ( £ {g* c ) = -2, the critical exponents fH£l[H3j 

-C«(ffc*) 10 -3d' ^ -Cfc*) 10 -3d 

Note that the presence of the dangerously irrelevant parameter l/a R precludes a direct 
calculation of the power laws precisely at the critical point (rather than approaching it 
from the active phase), and the derivation of 'hyperscaling' relations such as f3 = zua. 
Numerically, the PC critical exponents in one dimension have been determined to be 
v ~ 1.6, z m 1.75, a ~ 0.27, and (3 ~ 0.92 (HI E]- Perhaps not too surprisingly, 
the predictions (J123|) from the uncontrolled fixed- dimension expansion yield rather 
poor values at d = 1. Unfortunately, an extension to, say, higher loop order, is not 
straightforward, and an improved analytic treatment has hitherto not been achieved. 

6.5. BARW variants and higher-order processes 

• Levy flight BARW: 

Simulations clearly cannot access the PC borderline critical dimension d' c . This difficulty 
can be overcome by changing from ordinary diffusion to Levy flight propagation ~ D^p 17 . 
The existence of the power-law inactive phase is then controlled by the Levy exponent 
a, and in one dimension emerges for a > a c = 3/2 |SB] . 

• Multi-species generalizations of BARW: 

There is a straightforward generalization of the two-offspring BARW to a variant with 
q interacting species A4, according to A4 — > 3Ai (rate a), A4 A4 + 2Aj (j ^ i, 
rate a'), and A4 + A± —>■ only for particles of the same species. Through simple 
combinatorics aR_jo' R — > under renormalization, and the process with rate a' dominates 
asymptotically. The coarse-grained effective theory then merely contains the rate a' R , 
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corresponding formally to the limit q — > oo, and can be analyzed exactly. It displays 
merely a degenerate phase transition at a' c = 0, similar to the single-species even- 
offspring BARW for d > d' c , but with critical exponents v = 1/d, z = 2, a = d/2, and 
(3 = 1 = zva |148| Fl49j. The situation for q = 1 is thus qualitatively different from any 
multi-species generalization, and cannot be accessed, say, by means of a \ jq expansion. 

• Triplet and higher-order generalizations of BARW: 

Invoking similar arguments as above for k = 3, i.e., the triplet annihilation 3A — > 
coupled to branching processes, one would expect DP critical behavior at a phase 
transition with cr c > for any mmod3 = 1, 2. For m = 3, 6, . . ., however, special critical 
scenarios might emerge, but limited to mere logarithmic corrections, since d c = 1 in this 
case j!49| . Simulations, however, indicate that such higher-order BARW processes may 
display even richer phase diagrams jOj. 

• Fission / annihilation or the pair contact process with diffusion (PCPD) 

One may expect novel critical behavior for active to bsorbing state transitions if 
there is no first-order process present at all. This occurs if the branching reaction 
competing with A + A — > (0, A) is replaced with A + A — > (n + 2) A, termed 
fission/annihilation reactions in Ref. [2H], but now generally known as pair contact 
process with diffusion (PCPD) [21]. Without any restrictions on the local particle 
density, or, in the lattice version, on the site occupation numbers, the density obviously 
diverges in the active phase, whereas the inactive, absorbing state is governed by 
the power laws of the pair annihilation/coagulation process (231- By introducing site 
occupation restrictions, or alternatively, by adding triplet annihilation processes, the 
active state density becomes finite, and the phase transition continuous. In a field- 
theoretic representation, one must also take into account the infinitely many additional 
fission processes that are generated by fluctuations. Following Ref. |129j . one may 
construct the field theory action for the restricted model version, whence upon expanding 
the ensuing exponentials, see (|77jl. one arrives at a renormalizable action. Its RG 
analysis however leads to runaway RG trajectories, indicating that this action cannot 
represent the proper effective field theory for the PCPD critical point j2Hj. Since Monte 
Carlo simulation data for this process are governed by long crossover regimes, the 
identification and characterization of the PCPD universality class remains to date an 
intriguing open issue j21j. 

6.6. Boundaries 

In equilibrium critical phenomena it is well-known that, close to boundaries, the critical 
behavior can be different from that in the bulk (see Refs. |152[ 1153) for comprehensive 
reviews). As we will see, a similar situation holds in the case of nonequilibrium reaction- 
diffusion systems (see also the review in Ref. |154j ). Depending on the values of 
the boundary and bulk reaction terms, various types of boundary critical behavior 
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are possible. For example, if the boundary reaction terms ensure that the boundary, 
independent of the bulk, is active, while the bulk is critical, then we have the so-called 
extraordinary transition. Clearly, by varying the boundary/bulk reaction rates three 
other boundary transitions are possible: the ordinary transition (bulk critical, boundary 
inactive), the special transition (both bulk and boundary critical, a multicritical point), 
and the surface transition (boundary critical, bulk inactive). Defining r and r s as the 
deviations of the bulk and boundary from criticality, respectively, a schematic boundary 
phase diagram is shown in figure E3 In this review, for reasons of brevity, we will 
concentrate on the case of DP with a planar boundary |155t I156| I157j . Other cases 
(A + A — > with a boundary and boundary BARW) will be dealt with more briefly. 




Figure 10. Schematic mean field phase diagram for boundary DP. The transitions 
are labeled by 0=ordinary, E=extraordinary, S=surface, and Sp=special. 



• Boundary directed percolation: 

In this section, we will focus on the ordinary transition in boundary DP. The mean-field 
theory for this case was worked out in Ref . |157j , while the field theory was analyzed to 
one-loop order in Ref. |155j . 

As we have discussed earlier, the field theory for bulk DP is described by the action 
fjTBj) . Consider now the effect of a semi-infinite geometry {x = (xn,z), < z < oo}, 
bounded by a plane at z — 0. The complete action for bulk and boundary is then given 
by S = S c s + Sbd, where 

s» = f i*- 1 * (124) 

with the definitions s s = s(x\\,z = 0,t) and s s = s(x»,z = 0,t). This boundary term is 
the most relevant interaction consistent with the symmetries of the problem, and which 
respects the absorbing state criterion. Power counting indicates that the boundary 
coupling has scaling dimension [r s ] ~ k, and is therefore relevant. The presence of the 
wall at z = enforces the boundary condition 

d z s\ z=0 = r s s s . (125) 

This condition guarantees that a boundary term of the form sd z s is not required, even 
though it is marginal according to power counting arguments. 
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Since [r s ] ~ k the only possible fixed points of the renormalized coupling are 
A SR — > or ±00. Here we focus on the case r sR —>■ 00, corresponding to the ordinary 
transition. At this fixed point, the propagator in the presence of a boundary Go s can 
be written entirely in terms of the bulk propagator Go: 

G 0s (x\\,z,z',t) = G (x\\,z,z',t) - G (x\\,z,-z',t) . (126) 

Due to the above boundary condition, which implies that Go s (x\\, z, z', t)\ z= o = 0, we 
see that the appropriate boundary fields for the ordinary transition are not s s , s s , but 
rather s± = d z s\ z=0 , and s± = d z s\ z=0 . For example, in order to compute the order 
parameter exponent Pi at the boundary, defined by sr(z = 0,r) ~ \tr\^ {jr < 0), 
we must investigate how s± = d z s\ z= o scales. In mean- field theory, straightforward 
dimensional analysis yields f3\ = 3/2 [157] . Of course, to go beyond this simple mean- 
field picture and to incorporate fluctuations, we must now employ the machinery of the 
field-theoretic RG. 

Because of the presence of the surface, we expect to find new divergences which 
are entirely localized at the surface. These divergences must be absorbed into new 
renormalization constants, in addition to those necessary for renormalization of the 
bulk terms. At the ordinary transition, the new divergences can be absorbed by means 
of an additional surface field renormalization, yielding the renormalized fields: 

s ±R = Z Z l J 2 s ± , s ±R = Z Z l J 2 7s ± . (127) 

Note that the same factor Zq enters both renormalized surface fields, similar to the 
bulk field renormalization. The fact that one independent boundary renormalization is 
required translates into the existence of one independent boundary exponent, which we 
can take to be j3i, defined above. 

Consider now the connected renormalized correlation function G^' M ^ , composed of 
iV {s, s} fields and M {s±, fields. The renormalization group equation then reads 
(excepting the case iV = 0, M = 2 for which there is an additional renormalization): 

(4 - ^ <• - + ^ik + ^ + G ""' - ■ (l28) 

with the definitions (PT|) - ([Mj) and £ = nd K \nZ . Solving the above equation at the 
bulk fixed point using the method of characteristics, combined with dimensional analysis, 
yields 

Gp M \{x,t},D R7 T R} K,V* R ) ~ | Tfl |(iV + M), + MKl-,o) 6 (iV 1 M)^^ ) 0^^ ^ 

With e = 4 — d, and defining 7/0 = (o( v r) — e /12 + 0(e 2 ) (the value of the Co function 
at the bulk fixed point), we see that at the ordinary transition 

P 1 = P + v(l-Vo) = l-^ + 0(e 2 ) , (130) 

where we have used some results previously derived for bulk DP. The general trend 
of the fluctuation correction is consistent with the results of Monte Carlo simulations 
in two dimensions |156j and series expansions in d = 1 |158j . which give (3\ = 1.07(5) 
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and Pi = 0.73371(2), respectively. For dIP, an analogous analysis yields the boundary 
density exponent p 1 = 3/2 - lle/84 + 0(e 2 ) [12]. 

One unsolved mystery in boundary DP concerns the exponent T\ = zv — Pi, 
governing the mean cluster lifetime in the presence of a boundary |156j . This exponent 
has been conjectured to be equal to unity |159j : series expansions certainly yield a value 
very close to this (1.00014(2)) |158j . but there is as yet no explanation (field-theoretic 
or otherwise) as to why this exponent should assume this value. 

• Boundaries in other reaction-diffusion systems: 

Aside from DP, boundaries have been studied in several other reaction-diffusion systems. 
BARW (with an even number of offspring) with a boundary was analyzed using field- 
theoretic and numerical methods in Refs. |154 tll57(ll60j . As in the bulk case, the study 
of boundary BARW is complicated by the presence of a second critical dimension d' c 
which prevents the application of controlled perturbative e expansions down to d = 1. 
Nevertheless some progress could still be made in determining the boundary BARW 
phase diagram |157j . The situation is somewhat more complicated than in the case 
of DP, not only because the location of the bulk critical point is shifted away from 
zero branching rate (for d < d' c ), but also because the parity symmetry of the bulk 
can be broken but only at the boundary. The authors of Ref. |157j proposed that 
the one-dimensional phase diagram for BARW is rather different from that of mean- 
field theory: If a symmetry breaking A — > process is present on the boundary, 
then only an ordinary transition is accessible in d = 1; whereas if such a reaction is 
absent then only a special transition is possible. Furthermore, an exact calculation in 
d — 1 at a particular plane in parameter space allowed the authors of Ref. |157j to 
derive a relation between the Pi exponents at the ordinary and special transitions. It 
would be very interesting to understand this result from a field-theoretic perspective, 
but until controlled perturbative expansions down to d = 1 become possible, such an 
understanding will probably remain elusive. More details of these results can be found 
in Refs. fT54llT5Tj . 

Richardson and Kafri j!61l I162j analyzed the presence of a boundary in the simpler 
A + A — > reaction. For d < 2, they found a fluctuation- induced density excess develops 
at the boundary, and this excess extends into the system diffusively from the boundary. 
The (universal) ratio between the boundary and bulk densities was computed to first 
order in e = 2 — d. Since the only reaction occurring both on the boundary and in 
the bulk is the critical A + A —>■ process, this situation corresponds to the special 
transition. 

7. Open Problems and Future Directions 

As we have seen, enormous progress has been made over the last decade or so in 
understanding fluctuations in reaction-diffusion processes. Many systems are now rather 
well understood, thanks to a variety of complementary techniques, including mean-field 
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models, Smoluchowski approximations, exact solutions, Monte Carlo simulations, as 
well as the field-theoretic RG methods we have predominantly reviewed in this article. 
However, we again emphasize the particular importance of RG methods in providing 
the only proper understanding of universality. Despite these undoubted successes, we 
believe that there are still many intriguing open problems: 

• Already for the simple two-species pair annihilation process A + B — ► 0, field- 
theoretic RG methods have not as yet been able to properly analyze the asymptotic 
properties in dimensions d < 2 in the case of equal initial densities [2*T] . Moreover, the 
standard bosonic field theory representation appears not to capture the particle species 
segregation in multi-species generalizations adequately [108J. A viable description of 
topological constraints in one dimension, such as induced by hard-core interactions that 
prevent particles passing by each other, within field theory remains a challenge. 

• Branching-annihilating random walks (BARW) with an even number of offspring 
particles is still poorly understood in d = 1, due to the existence of the second critical 
dimension d' c ,148] 1149] . A systematic extension of the one-loop analysis at fixed 
dimension to higher orders has not been successfully carried out yet. Ideally one would 
like to find a way of circumventing this difficulty, in particular to understand why certain 
one-loop results (for the exponent (3 and the value of d' c |85j) appear to be exact, even 
when the two-loop corrections are known to be nonzero. Nonperturbative numerical 
RG methods might be of considerable value here j!5()[ 1151] . There is also an interesting 
suggestion for a combined Langevin description of both DP and PC universality classes 

163 , but the ensuing field theory has yet to be studied by means of the RG. 

• Despite intensive work over recent years the status of the the pair contact process with 
diffusion (PCPD) 23j is still extremely unclear. In particular, even such basic questions 
as the universality class of the transition, remain highly controversial. Since simulations 
in this model have proved to be very difficult, due to extremely long crossover times, 
it appears that only a significant theoretical advance will settle the issue. However, 
the derivation of an appropriate effective field theory remains an unsolved and highly 
nontrivial task [25J. Other higher-order processes also appear to display richer behavior 
than perhaps naively expected 0. 

• Generally, the full classification of scale-invariant behavior in diffusion-limited 
reactions remains a formidable program, especially in multi-species systems; see 
Refs. jHJ E| for an overview of the current data from computer simulations. To date, 
really only the many-species generalizations of the pair annihilation reaction as well as 
the DP and dIP processes are satisfactorily understood. 

• An important, yet hardly studied and less resolved issue is the effect of disorder in 
the reaction rates, especially for active to absorbing state transitions. A straightforward 
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analysis of DP with random threshold yields runaway RG flows |164| . which seem to 
indicate that the presence of disorder does not merely change the value of the critical 
exponents, but may lead to entirely different physics (see, e.g., Ref. (165 ). This may 
in turn require the further development of novel tools, e.g., real-space RG treatments 
directly aimed at the strong disorder regime |166^ I167j . 

• In contrast with the many theoretical and computational successes, the subject of 
fluctuations in reaction-diffusion systems is badly in need of experimental contact. Up 
to this point, the impact of the field on actual laboratory (as opposed to computer) 
experiments has been very limited. In this context, the example of Directed Percolation 
(DP) seems especially relevant. DP has been found to be ubiquitous in theory and 
simulation, but is still mostly unobserved in experiments, despite some effort. Ideally, 
one would like to understand why this is the case: could it be due to disorder or to the 
absence of a true absorbing state? 

• There are a number of additional extensions of the field-theoretic approach presented 
here that could further improve our understanding of reaction-diffusion systems. For 
example, Dickman and Vidigal have shown how to use this formalism to obtain the 
full generating function for the probability distribution of simple processes [168 ; Elgart 
and Kamenev have used the field theory mapping to investigate rare event statistics 

169 ; and Kamenev has pointed out its relation to the Keldysh formalism for quantum 
nonequilibrium systems |17()j . Path-integral representations of stochastic reaction- 
diffusion processes are now making their way into the mathematical biology literature 
[T7THT72] . 

We believe that these questions and others will remain the object of active and fruitful 
research in the years ahead, and that the continued development of field-theoretic RG 
methods will have an important role to play. 
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